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ARRAY  SIGNAL  PROCESSING  UNDER  MODEL  ERRORS  WITH 
APPLICATION  TO  SPEECH  SEPARATION 


Abstract 

by 

Victor  C.  Soon 


In  recent  years  there  have  been  increasing  interest  in  eigenstructure  based  array 
signal  processing  techniques  especially  in  the  areas  of  radar  and  sonar  processing 
and  in  signal  separation  applications  as  found  in  communications  or  speech  acquisi¬ 
tion/enhancement.  However,  model  errors  are  not  accounted  for  in  these  methods. 
Unfortunately,  the  problem  of  model  error  is  a  ubiquitous  part  of  real  systems;  they 
arise  from  such  phenomena  as  non-ideal  sensor  characteristics  (eg.,  sensor  calibration 
errors  in  phase  and  gain,  etc.),  unknown  source  characteristics  and  locations,  etc. 
This  report  is  broadly  divided  into  two  parts.  The  first  part  examines  the  effects  of 
model  errors  on  the  performance  of  existing  eigenstructure  based  methods  in  array 
signal  processing.  A  statistical  analysis  of  the  ESPRIT  (Estimation  of  Signal  Param¬ 
eters  via  Invariance  Techniques)  algorithm  under  random  model  errors  is  performed. 
The  analysis  provides  interesting  insight  into  the  sensitivity  of  ESPRIT  to  model 
errors;  in  particular,  for  uniform  linear  arrays  of  sensors,  it  is  found  that  the  mean 
square  error  of  DOA  (direction  of  arrival)  estimates  found  using  ESPRIT  is  almost 
totally  dependent  on  errors  in  sensor  phases  and  not  that  of  sensor  gains.  The  second 
part  of  the  report  deals  with  the  development  of  algorithms  that  take  model  errors 
into  account.  Two  approaches  are  proposed.  The  first  is  based  on  a  signal  subspace 
constraint  on  the  possible  set  of  model  parameters  as  determined  from  the  array  co- 
variance.  This  constraint  is  incorporated  into  an  iterative  procedure  which  calibrates 


the  axray  under  model  errors  and  unknown  source  signals.  The  second  approach 
incorporates  blind  identification  and  clustering  for  the  source  estimation  problem  un¬ 
der  model  errors  or  uncertainties.  The  proposed  approach  is  shown  to  be  robust  to 
the  silmutaneous  presence  of  such  uncertmnties  such  as  unknown  sensor  and  chamnel 
gains,  unknown  combinations  of  near-field  and  far-field  sources,  unknown  combina¬ 
tions  of  narrowband  and  wideband  sources,  unknown  source  spectral  characteristics 
and  unknown  number  of  sources. 
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INTRODUCTION 

TP he  processing  of  signals  obtained  from  arrays  of  sensors  is  of  interest  in  a  wide 
range  of  problems  in  radar  and  sonar  processing  where  the  presence  or  absence  of 
source  signals,  their  direction-of-arrivals  (DOA),  signal  power,  etc.  axe  of  interest  (see 
eg.  [53,  50,  83],  in  speech  recognition  and  teleconference  applications  where  source 
signals  from  spatially  distributed  locations  are  to  be  separated  from  mixed  measure¬ 
ments  (the  so  called  “cocktail  party  problem”),  see  eg.  [25,  36,  32,  12,  16,  57,  22).  In 
communications  one  may  be  interested  in  extracting  a  desired  signed  from  jammer- 
contaminated  measurements,  see  eg.,  [11].  In  biomedical  signal  processing  a  set  of 
measurements  taken  from  an  array  of  electrodes  may  be  processed  to  extract  signals 
from  individual  neuronal  firings,  etc.  The  rapidly  expanding  demand  for  mobile  com¬ 
munications  such  as  in  cellular  telephone  technology  provides  an  interesting  area  of 
possible  application  for  array  signal  processing  techniques  where  their  capability  of 
resolving  signals  based  on  their  spatial  locations  may  make  it  possible  to  increase 
channel  capacity,  see  eg.,  [66].  In  radio  astronomy,  arrays  of  radio  observatories  are 
used  in  obtaining  higher  resolution  maps  of  galaxies  and  quasars.  Radar  imaging  is 
also  one  application  where  array  signal  processing  is  used. 

In  most  of  these  applications  a  common  model  is  assumed  to  adequately  model  the 
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physical  processes  in  the  ipplications.  When  assumptions  of  the  techniques  in  array 
signal  processing  are  satisfied,  these  techniques  are  capable  of  providing  superior  per¬ 
formance  as  compared  to  conventional  methods.  This  is  the  case  in  DOA  estimation 
where  the  modern  eigenstructure  based  methods  are  capable  of  resolving  DOA  an¬ 
gles  separations  that  are  much  smaller  than  the  conventional  (classical)  beamforming 
method  which  is  limited  to  the  Rayleigh  resolution  limit.  These  techniques  however 
are  dependent  on  model  assumptions  that  may  not  be  true  in  actual  systems,  i.e., 
model  errors  are  not  accounted  for  in  these  methods.  Unfortunately,  the  problem  of 
model  error  is  a  ubiquitous  paurt  of  real  systems:  they  arise  from  such  phenomena 
as  non-ideal  sensor  characteristics  (eg.,  sensor  calibration  errors  in  phase  and  gain, 
etc,),  errors  in  the  locations  and  orientation  of  the  sensors  within  the  array,  source 
characteristics  that  deviate  from  the  assumed  source  model  (eg.,  sources  are  not  ideal 
point  sources,  intensity  of  sources  may  favor  certain  direction  of  propagation  resulting 
in  unequal  intensities  observed  at  different  sensors),  channel  ii?homogenities,  sources 
in  the  near-held,  etc. 

The  goal  of  this  report  is  thus  to  study  the  problem  of  array  signal  processing 
under  model  uncertainties.  A  performance  analysis  of  the  ESPRIT  (Estimation  of 
Signal  Parameters  via  Rotational  Invariance  Techniques)  ,  a  popular  eigenstructure 
based  method,  under  random  model  uncertainties  is  undertaken  to  further  understand 
its  robustness  to  model  errors.  Two  techniques  are  then  proposed  for  this  problem: 
the  first  is  based  on  a  signal  subspace  constraint  on  the  possible  sets  of  model  errors 
and  DOA  angles  while  the  second  is  based  on  the  incorporation  of  blind  identification 
techniques  with  clustering  techniques.  The  integration  of  these  techniques  and  ideas 
may  make  it  possible  to  apply  array  signal  processing  under  less  ideal  conditions. 

Chapter  2  presents  the  model  formulation  for  the  array  signal  processing  problem 
and  its  areas  of  application.  Chapter  3  covers  the  basic  mathematical  machinery 
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needed  for  the  study  in  the  report.  Chapter  4  gives  a  broad  overview  of  various 
aspects  of  array  signal  processing,  from  conventional  methods  to  eigenstructure  based 
methods  and  other  aspects  such  as  model  order  estimation  and  extensions  fot  coherent 
and  wideband  situations.  Chapter  5  presents  the  performance  analysis  of  ESPRIT 
under  random  model  errors.  Chapter  6  contains  the  signal  subspace  approach  to  array 
signal  processing  under  model  errors.  Chapter  7  presents  an  approach  bcLsed  on  the 
incorporation  of  blind  identification  and  clustering  techniques  to  solve  the  problem. 


» 
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MODEL  FORMULATION 

In  this  chapter,  we  shall  cover  the  basic  model  formulation  of  the  array  signal  pro¬ 
cessing  problem.  .Narrowband  and  wideband  formulations  are  then  discussed.  Model 
parameters  such  as  time  delays,  direction-of-arrival  angles  of  source  signals  and  chan¬ 
nel  coefficients  are  then  looked  at.  In  the  final  section  areas  of  application  for  array 
signed  processing  methods  are  examined. 

2.1  Notations  and  Preliminaries 

Throughout  this  thesis,  we  assume  the  following  notational  convention:  All  matrix 
quantities  w’ill  be  denoted  by  upper  case  bold  faced  symbols  while  vector  quantites 
will  be  denoted  by  lower  case  bold  faced  symbols.  Also, 

•  {.)‘  -  Transpose  of  (.) 

•  {.)^  -  Complex  conjugate  transpose  of  (.) 

•  {.)”*  -  Inverse  of  (.) 

•  (.)^  -  Pseudo- in  verse  of  (.) 

•  ii  M  IIf  —  Frobenius  norm  of  M 
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•  (I  V  |j  —  Euclidean  norm  of  the  vector  v. 


•  II  M  II2  -  2- norm  of  matrix  M 

•  m  —  Number  of  sensors  within  array 

•  n  —  Number  of  sources 

•  —  Location  of  i-th  sensor  in  array 

•  Tifc  —  Time  delay  of  k-th.  source  signal  at  i-th  sensor. 

•  6k  —  Direction-of-arrival  of  k-th  source  signal  measured  relative  to  y-axis. 

VVe  write  the  general  formulation  of  the  problem  in  the  following  manner.  Let  the 
i-th  sensor  output  be  denoted  by 

General  Model 

n 

Xi{t) -'^aikSkit i  =  l,2,---,m.  (2.1) 

where 

Cik  :  channel  coefficient 

Sk  :  L'-th  source  signal 

T}i{t)  ;  additive  noise  at  i-th  sensor. 

Tik  :  Time  delay  at  i-the  sensor  of  fc-th  source. 

Fig.  2.1  shows  the  model. 

The  channel  coefficients  are  the  result  of  the  sensor  and  source  characteristics  - 
i.e.  the  sensors  and  sources  directionality,  sensitivity,  gains,  phases,  etc.  and  channel 
parameters,  inhomogeneities,  etc.  The  time  delays  are  due  to  the  spatial  distribu¬ 
tion  of  the  sensors  and  locations  of  the  source  signals  -  they  are  functions  of  sensor 
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Figure  2.1:  Array  Signal  Processing  Model 

locations,  source  locations,  distance  between  sensors,  etc.  Propagation  of  the  sign2ds 
through  the  channel  can  be  far-field  or  near  field  or  any  combination  thereof. 

The  general  problem  which  is  wideband  in  form  is  such  that  the  channel  and  signal 
components  cannot  be  separated  in  the  time-domain.  However,  under  narrowband 
conditions  (which  will  be  discussed  in  the  next  section),  Sk{t)  may  be  cissumed  to  be 
baseband  and  the  model  written  as: 

Narrowband  Model 

x(t)  =  As(t) -j- ri(t)  (2.2) 

■  ‘ 

A  =  .  .  (2.3) 

a  -—ju*oTmi  ...  fj 

where  ujq  is  the  center  frequency  of  the  signals. 
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The  specific  form  of  the  parameterization  of  the  time  delays  will  be  elaborated  on 
in  the  sequel. 

The  following  assumptions  will  be  assumed  to  hold  for  all  the  coming  discussions 
in  this  thesis. 

Assumptions 

(Al)  The  matrix  A  in  Eqn.  (2.3)  has  full  column  rank. 

(A2)  s(t)  is  zero  mean  stationary  with  a  positive  definite  covariance  matrix. 

(A3)  Tj(t)  is  additive  white  Gaussian  noise,  mutually  uncorrelated  with  s(t). 

(A4)  The  number  of  sources,  n  is  assumed  knowm  or  estimated  previously. 

2.2  Narrowband  and  Wideband  Formulations 

As  from  Eqn.  (2.1) 

n 

H  -  T.fc),  t  =  1, 2,  •  •  • ,  m.  (2.4) 

fe=i 

If  sjt(t)  is  bandlimited  with  center  frequency  wo,  we  write 

-Sjt(0  =  Sklt)cos(uot  +  <^k(i)) 

where  Sfe(t)  denotes  the  envelope  of  Sk(t)  and  <j>k{t)  denotes  the  phase  associated  with 

Sk{t). 

Under  narrowband  assumptions  on  Sk{t),  the  time  delays  r,>  are  small  enough 
such  that  Sfc(f)  ~  Sk{t  -  r.fc)  and  <i>k{t)  ~  (f>k{t  -  t,*). 

Definition:  A  baseband  signal  Skit)  is  said  to  be  narrowband  if  the  reciprocal  of 
its  highest  frequency  component  is  much  larger  than  the  maximum  delay  time  of  the 
signal  propagating  through  the  sensor  array  in  question,  see  also  [78]. 
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The  signal  is  said  to  be  wideband  if  it  is  not  narrowband. 

Thus  by  demodulating  the  received  signal  at  the  i  —  th  sensor  down  to  the  baseband 
level  through  a  quadrature  demodulator  and  passing  the  results  through  low  pass 
filters,  we  can  write  the  outputs  as: 

n 

In  —  Phase  —  Component  :  Xi{t)  =  ^  a,fccos(a;or,jt  +  <i>ic{t))sk{t)  (2.5) 

k=l 

n 

Quadrature  —  Phase  —  Component :  i,(t)  =  —  aiksin{ujQTik  +  <i>k{t))sk{t)  (2.6) 

fc=i 

Combining  both  the  in-phase  and  quadrature-phase  components  we  can  write  the 
complex  valued  i-th  sensor  output  as 

1.(0  =  -j2  a.;ke-^‘*^"-*i*(0e-‘^‘‘“  (2.7) 

fc=t 

Remark:  The  result  of  the  quadrature  de-modulation  of  a  real  valued  white  gaussian 
noise  is  circular  complex  valued  uncorrelated  gaussian  noise,  see  [35]. 

Writing  the  array  observations  as  a  vector,  we  get  Eqn.(2.2) 

x(0  =  As(0  +  »7(0  (2-8) 

Here  s(0  =  [si(0e'”'''^'^‘*  •  •  •  ■s„(0e~-''*’"^‘^]^  is  the  source  signal  vector  and  77(0  = 
[771(0  ■  •  •  ^m(0]  ^he  additive  observation  noise  vector.  Also, 

••• 

where  wq  is  the  center  frequency  of  the  signals.  Fig.2.2  shows  a  block  figure 
diagram  of  the  quadrature  demodulator. 

Under  wideband  conditions,  the  assumption  sjt(0  —  ^k{t  —  r^k)  and  ^k(t)  —  <i>k{t  — 
Tik)  are  no  longer  valid  and  the  array  observations  may  no  longer  be  written  in  the 


COS(0[^t) 


►In-Phtse 


>Quad-Phue 


Figure  2.2;  Quadrature  Demodulator  Block  Diagram 

form  of  Eqn,  (2.8),  at  least  not  in  the  time-domain.  However,  by  taking  the  Fourier 
Transform  of  the  general  array  observations,  we  get 

x(w)  =  A(u;)s(u>)  +  ij{u)  (2.10) 


•  •  • 

••• 

(2.11) 


where  u  denotes  frequency  in  radians. 


The  matrix  A  is  usually  referred  to  as  the 


steering  matrix  in  the  literature. 


2.3  Time  delays,  Tik  and  Direction-of- arrivals  (DOA) 
Ok  Relations 


In  this  section,  the  form  of  the  time  delays  Tik  will  be  discussed.  Due  to  the  spatial 
distributions  of  the  sources  and  sensors,  different  sensors  will  observe  the  same  signal 
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at  different  time  delays.  Thus,  by  selecting  the  location  of  the  reference  sensor  as 
(0,0),  the  A:-th  source  signal  will  be  observed  at  the  time  delay  value  r,*  at  the  i- 
th  sensor.  The  direction  of  arrival  of  the  k-ih  source,  denoted  by  6k,  is  defined  with 
reference  to  the  vertical  (u)  axis  at  the  reference  sensor.  Without  any  loss  of  generality 
one  can  select  the  reference  sensor  location  as  the  origin. 

There  are  two  cases  of  interest. 

Far«Field:  Consider  the  far-field  case,  i.e.,  the  plane  wave  propagation  assump¬ 
tion  is  valid,  see  Fig.  2.3.  From  Fig.2.3,  it  can  be  easily  shown  through  simple 
analytical  geometry  that. 


Tik  =  cos  Ok  +  fli  sin  } 
c 


where  c  denotes  the  velocity  of  propagation  of  the  signal. 


(2.12) 


Example  (Uniform  Linear  Array):  The  uniform  linear  array  (ULA)  is  an 
important  special  case  in  array  signal  processing.  It  corresponds  to  the  case  where 
ther  sensors  are  aligned  in  a  single  axis  and  spaced  with  equal  distance  d  between 
adjacent  sensors,  see  Fig.2.4.  In  this  case,  the  time  delays  may  be  written  as: 


{i  —  l)dsin  6k 


(2.13) 


When  the  channel  coefficients  a,>  are  all  equal  to  1.0,  the  form  of  the  steering  matrix 
A  from  (2.11)  is  written  as 


g-JU/T„(0) 

g-Ju;T„(l) 


(2.14) 


1)  g— JwTB(m-l) 


where  r*  =  isBh.  T^e  steering  matrix  A  in  this  case  has  a  special  form  called 
Vandermonde.  It  has  been  shown  that  it  has  full  column  rank  when  the  {r*}  terms 


are  distinct. 
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Figure  2.3:  Far-Field  Propogation 

Near-Field:  For  the  case  where  the  far-field  assumption  does  not  hold  the 

wavefronts  cannot  be  assumed  to  be  planar  but  rather  are  spherical.  Fig.  2.5  shows 
the  situation  in  the  near  field.  The  form  for  the  time  delay  is  then  written  as 

Tik  =  ^.^(far  -  field)  -  +  i',  cos^k)^  - 

where  T,k(far  —  field)  is  the  time  delay  assuming  far-field  propagation  and  R  is  the 
range  or  distance  of  the  source  signal  from  the  reference  sensor,  see  also  [17].  From 
Eqn.(2.15)  it  can  be  seen  that  the  near-field  time  delay  has  a  correction  term  (1//2) 
due  to  the  curvature  of  the  propagating  wave.  It  is  obvious  that  as  the  range  gets 
very  large,  the  time  delay  approaches  the  far-field  form. 
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(0.0) 

(d.O) 

(2d.0) 

((m-l)d,0) 

Sensor  Amy 

Figure  2.4:  Uniform  Linear  Array  (ULA) 

2.4  Channel  Coefficients  Relations 

The  channel  coefficients  a,fc  terms  are  a  result  of  the  combined  effects  of  the  chan¬ 
nel  attenuation  characteristics,  the  sensor  characteristics  gains  and  the  geometrical 
attenuation  due  to  wave  propagation. 

Channel  Attenuation:  The  channel  attenuation  characteristics  are  dependent 
on  the  type  of  channel  in  question,  for  example  it  may  be  the  atmospheric  chan¬ 
nel  for  radar,  underwater  sea  channel  for  sonar  or  the  air  sound  channel  for  speech 
applications.  The  physical  process  by  which  a  source  signal  is  attenuated  through 
these  channels  are  often  highly  complex  and  difficult  to  model.  In  underwater  sonar 
channels,  signal  attenuation  is  caused  by  a  combination  of  molecular  relaxation  pro¬ 
cesses  (such  as  chemical  relaxation  in  salt  water)  and  frictional  heat  losses  due  to  the 
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Source 


\7 

Figure  2.5:  Near-Field  Propagation 

relative  motion  of  portions  of  the  transmission  medium,  see  [9,  49]. 

Sensor  Gains:  The  sensor  gains  are  due  to  the  combined  effects  of  the  sensor 
receiving  characteristics,  the  gains  of  the  amplifiers  and  the  subsequent  conditioning 
done  on  the  signals  such  as  quadrature  de-modulating,  etc. 

Geometrical  Attenuation:  The  physics  of  wave  propagation  in  electromagnetic 
and  acoustic  theory  is  governed  by  the  wave  equation  (a  d’Alembertian  equation), 
with  vector  quantities  in  electromagnetism  and  scalar  quantities  in  acoustics.  When 
a  spherical  propagation  mode  is  assumed  in  the  wave  equation,  the  intensity  of  the 
fc-th  signal  at  the  f-th  sensor  is  inversely  proportional  to  Rik  which  is  the  radial 
distance  of  the  k-ih.  source  to  the  i-th  sensor,  see  for  example  [55,  49).  Thus  the 
channel  coefficient  are  inversely  proportional  to  Rik.  Generally,  it  is  clear  that  the 
channel  coefficients  are  not  equal  to  each  other  since  Rik  s  are  not  the  same.  However, 
in  the  case  where  the  dimensions  of  the  array  are  small  in  comparison  to  the  distances 
Rik,  the  geometrical  attenuations  may  be  assumed  to  be  approximately  equal. 
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2.5  Areas  of  Interest  and  Application  of  Model 


The  general  model  and  the  narrowband  model  of  Eqns.(2.10},(2.8)  and  the  associated 
parameters  are  a  focus  of  interest  in  many  applications.  For  example,  the  time  delays 
and  DOAs  of  the  various  source  signals  are  of  interest  in  radar  and  sonar  processing 
while  the  source  signals  themselves  are  of  interest  in  communications  applications. 

2.5.1  Array  Signal  Processing  for  DOA  Estimation 

In  this  application,  the  DOA  of  various  source  signals  are  of  interest.  When  the 
channel  coefficients  are  known  ,  or  normalized  to  unity,  or  the  sensors  in  the  array 
are  identical  the  models  can  be  further  simplified  to  the  following. 

Wideband  Model 

x(w)  =  A(w)s(u;)  +  T}{u)  (2.16) 


A  = 


g-JWTSl 


where  u  denotes  frequency  in  radians. 


Narrowband  Model 


g—JU/Tfnn 


x(0  =  As(<)  +  T]{t) 


(2.17) 


(2.18) 


g-JWDTll 

g-JW0Tl„ 

g-ju«oTj| 

g— juJoTjn 

,  , 

g*Ju^Tn,n 

(2.19) 


where  ujq  is  the  center  frequency  of  the  signals. 
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Thus  given  the  observation  of  x(t)  over  N  snapshots  or  samples,  the  objective  is 
to  estimate  the  DOA  angles  for  k  = 


2.5.2  Spectrum  Estimation 


There  is  a  duality  between  the  estimation  of  multiple  complex  sinusoids  in  noise  and 
the  estimation  of  DOAs  using  uniform  linear  arrays.  The  scalar  observation  can  be 
written  as 

=  +  >,(()  (2.20) 

i=l 

where  w,  denotes  the  frequency  (in  radians)  of  the  i-th  sinusoid,  <f>i  denotes  the  phase 
of  the  i-th  sinusoid  and  s,  denotes  the  intensity  of  the  i-th  sinusoid.  Assuming  that 
the  observation  is  sampled  uniformly  at  the  sampling  frequency  of  /,,  the  discrete 
scalar  observation  is  written  as 

x{k)  =  ^  ^  (2.21) 

«ssl 

By  collecting  m  consecutive  discrete  observations  into  a  column  vector  form  we  get 

x{k)  =  As{k)  +  r){k)  (2.22) 


where  s  =  S2e-'^  •  •  •  T){k)  =  [T}{k)  r}{k  —  1)  ■  •  •  r){k  —  m  -f  1)]‘  and 


g-ju/»(o)//,  g-jwnio)//,  • 

... 

•  •  •  •  • 

g-jwt(m-t)/f,  ...  g->w„(m-l)//. 


(2.23) 


Thus  the  matrix  A  is  in  the  same  form  as  in  the  uniform  linear  array  case  discussed 


in  a  previous  section,  i.e.,  it  is  Vandermonde. 


2.5.3  Signal  Separation  Applications 

The  extraction  of  certain  desired  source  signals  from  contaminated  or  jammed  obser¬ 
vations  can  be  found  in  areas  such  as  speech  or  tele-conferencing  processing  where 


individual  speech  signals  must  be  extracted  from  an  environment  with  background 
and  other  noise  sources,  in  communications  where  sources  of  interference  may  degrade 
the  received  signals  and  in  EMG  medical  processing  where  the  evoked  potentials  from 
various  neurons  must  be  separated.  The  model  may  be  either  narrowband  or  wide¬ 
band  depending  on  the  specific  application  of  interest.  Certainly  in  communications 
systems  which  employ  narrowband  signals,  the  narrowband  model  is  appropriate.  In 
speech  applications  however,  the  wideband  model  is  a  better  representation  of  the 
problem  since  speech  signals  are  wideband  in  nature. 

Remark:  Note  that  the  emphasis  in  signal  applications  is  on  the  extraction  of  indi¬ 
vidual  signals  from  observations  where  they  are  mixed.  This  requires  the  estimation 
of  the  channel  coefficients  and  time  delays  {r,*.}. 

2.5.4  Echo  Resolution  of  Signals  with  Known  Shapes 

The  problem  of  resolving  echoes  of  multiple  signals  of  known  shape  may  also  be 
framed  to  fit  the  model  of  (2.18).  This  is  achieved  by  a  suitable  parameterization  of 
the  matrix  A  in  terms  of  the  known  shapes,  see  [7].  Suppose  that  a  signal  g{t)  of 
known  shape  is  used  to  probe  a  medium  with  n  scatterers.  n  echoes  of  this  signal 
would  be  reflected  back  to  the  receiving  sensor  and  we  write  the  received  sensor  output 
as 

=  +  (2.24) 

1=1 

where  {s,}  models  the  random  gains  due  to  channel  propagation  and  scatterer  char¬ 
acteristics  and  r,  denotes  the  time  delay  of  the  i-th  echo.  Next,  suppose  that  the 
signal  g{t)  is  sent  to  probe  the  channel  repetitively,  say  N  times  and  assume  that  the 
echo  returns  would  have  died  out  at  the  beginning  of  each  probe  period.  Then  there 
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is  an  ensemble  of  N  sensor  observation  periods,  which  is  written  as 

|x{0  =  +  a-  =  l,2,  --,.V  and  t  =  (2.25) 

where  Si{k)  denotes  the  i-th  echo’s  random  gain  during  the  k-th.  observation  interval 
and  t  denotes  the  discrete  time  index. 

By  re-writing  the  scalar  observation  into  a  vector  form,  we  get 

x{k)  =  As(it)  -h  r}(k)  (2.26) 

where  x  =  [x(l)  x(2)  •  •  •  x(m)]‘,  s  =  [si  sq  •  •  ■  v  =  [t/(1)  tj(2)  •  •  •  t/(m))‘.  Here 
define  the  ’steering  matrix'  A  =  -  n)  1  g(<  -  ^2)  1  ■  •  •  1  g(<  -  Tn)],  with  g(f)  = 

[p(l)  5(2)  Thus  the  model  in  this  case  fits  the  narrowband  model  with 

the  steering  matrix  A  parameterized  in  terms  of  the  signal  g{t)  of  known  shape. 
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MATHEMATICAL 

PRELIMINARIES 

In  this  chapter  the  basic  mathematical  results  of  matrix  and  perturbation  theory 
is  presented.  Two  matrix  decompositions  are  covered:  the  eigenvalue  eigenvector 
decomposition  (EED)  and  the  singular  value  decomposition  (SVD).  Orthogonal  pro¬ 
jectors  and  solutions  to  the  deterministic  least  square  problem  are  presented.  The 
chapter  ends  with  basic  perturbation  results  o.i  the  perturbations  of  eigenvalues  and 
eigenvectors  of  a  Hermitian  covariance  matrix  and  the  perturbation  of  a  pseudoin¬ 
verse. 


3.1  Matrix  and  Vector  Norms 

Throughout  this  thesis  we  will  utilize  the  following  aorms  for  vectors  and  matrices. 

Definitions:  The  Euc'idean  norm  of  a  column  vector  v  of  dimension  m  is  defined 
as 

m 

II  V  11=  (E  I  I’)’"  (3-1) 

«=l 

Thus  the  Euclidean  norm  is  nothing  but  the  length  of  the  vector  v. 

The  2-norm  of  the  matrix  M  is  defined  as  that  induced  by  the  Euclidean  vector 
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nonii  as  in  the  following 

II  M  ||2=  sup  I!  Mv  II  (3.2) 

where  the  supremum  is  taken  over  all  vectors  v  of  unit  norm. 

Finally,  the  Frobenius  norm  of  a  matrix  M  =  [mi*]  is  defined  as 

II M  ilf=  (E  E  I  "-i*  (3-3) 

i=l  fe=l 

The  norms  defined  above  are  also  consistent,  i.e., 

||MN11<1|M|1||N||  (3.4) 

I|Mv||<1|M|||1  v||  (3.5) 

Next  we  state  certain  observations  regarding  the  matrix  theory. 

3.2  Eigenstructure  Decompositions 

In  array  signal  processing,  the  decomposition  of  a  covariance  matrix  is  frequently  of 
interest.  The  covariance  matrix  is  defined  as  R  =  E{xx^}.  Thus  it  is  Hermitian  and 
positive  semi-definite. 

Given  a  Hermitian  matrix  R,  there  are  several  different  decompositions  that  are 
available.  In  this  thesis,  we  will  be  concerned  with  the  Eigenvalue- Eigenvector  De- 
compostion  (EED)  and  the  Singular  Value  Decomposition  (SVD),  see  [15]. 

Eigenvalue-Eigenvector  Decomposition:  Given  a  Hermitian  matrix  R  of  size 
m-by-m  and  rank  n,  where  m  >  n,  there  exists  a  decomposition  of  R  in  the  following 
manner 


R  =  UAU^ 


(3.6) 


where  U  is  a  unitary  matrix  (i.e..  UU^  =  U^U  =  I)  and  A  is  a  diagonal  matrix 
consisting  of  the  real-valued  eigenvalues  of  R,  i.e., 

■A,  O' 

A=  (3.7) 

.0  A,„  . 

Example:  Consider  the  covariance  matrix  of  the  observation  vector  as  found  in  the 
narrowband  model  of  Eq.(2.8)  and  apply  assumptions  A1-A4  of  the  previous  chapter. 
Thus 

R  =  E{xx^}  =  AQA^  -f  (3.8) 

where  A  has  full  column  rank  (from  assumption  Al),  Q  is  the  positive  definite  co- 
variance  matrix  of  the  source  signals  (from  assumption  A2)  and  cr^  is  the  variance  of 
the  additive  white  Gaussian  noise  vector  (from  assumption  A3). 

Its  eigenvalue-eigenvector  decomposition  can  then  be  found  as 

R  =  UsAsUs^  +  <T*UnUn^  (3.9) 

where  U  =  (U,U„]  is  unitary,  As  is  an  n-by-n  diagonal  matrix  of  eigenvalues  greater 
than  <7^.  The  subspaces  spanned  by  the  column  vectors  of  U,,  U„  are  referred  to  as 
signal  and  null  (or  noise)  subspace,  respectively. 

Similarly,  in  the  case  of  the  wideband  model  of  Eq.(2.10) 

R(u;)  =  £;{x(u;)x^(u;)} 

=  A(c4;)Q(u;)A^(a;)  -|-  cr^(u;)I 

=  Us(u;)As(u;)U*^(u;)  a^(u;)U„(w)UnH^)  (3.10) 

Singular  Value  Decomposition:  Given  a  rectangular  matrix  A  of  size  m-by-n 
2Uid  rank  n,  where  m  >  n,  the  Singular  Value  Decomposition  of  A  is: 


A  =  UEVt 


(3.11) 


where  U.  V  are  unitary  matrices  of  sizes  m-by-m  and  n-by-n,  respectively  and  S  is 


S  = 


where  cr,,  i  =  1, 2, n  are  real. 


(3.12) 


Remark:  The  SVD  and  EED  of  a  positive  definite  Hermitian  matrix  (i.e.,  con¬ 
jugate  symmetric)  matrix  are  the  same. 


3.3  Orthogonal  Projectors 

Definition:  A  square  matrix  P  is  said  to  be  an  orthogonal  projectoronto  the  subspace 
S  if  it  satisfies  the  following: 

(I)  =  P,  it  is  idempotent 

(II)  Range{P)  =  S 

(III)  P^  =  P,  it  is  Hermitian. 

Remarks: 

(i)  An  orthogonal  projector  onto  a  given  subspace  is  unique,  see  (15|. 

(ii)  The  orthogonal  projector  P-*-  which  projects  onto  the  subspace  orthogonal  to 
that  of  P  satisfies  P-*-  -f-  P  =  I. 

Example:  Consider  the  covariance  matrix  R  =  P{xx^}  =  AQA^  +  <t^I  as  in 
Eq.(3.8).  Its  eigenvalue-eigenvector  decomposition  is  then 

R  =  UsAsUs^  -f  (T^UnUn^  (3.13) 

where  U  =  [U,Un]  is  unitary,  and  Ag  is  m  n-by-n  diagonal  matrix  of  eigenvalues 
greater  than  <7^.  The  subspaces  spanned  by  the  column  vectors  of  U,,  U„  are  referred 
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to  as  signal  and  null  (or  noise)  subspace,  respectively.  The  orthogonal  projector  for 
the  signal  subspace  may  be  defined  as: 

P,  =  (3.14) 

The  above  projector  may  be  shown  to  satisfy  conditions  (I-III)  above.  Similarly  the 
orthogonal  projector  for  the  noise  subspace  may  be  shown  to  be 

P«  =  UnUn^  (3.15) 

Since  the  eigenvector  matrix  U  is  unitary,  it  follows  that  the  signal  subspace  orthog¬ 
onal  projector  is  also  orthogonal  to  the  noise  subspace  orthogonal  projector,  (i.e.. 
P,  ±  P„)  and  that 

P,+Pn  =  I  (3.16) 

Example:  Consider  a  matrix  A  of  full  column  rank  and  dimensions  m-by-n.  The 
orthogonal  projector  P^  for  the  subspace  spanned  by  the  column  vectors  of  A  can 
be  shown  to  be 

P^  =  A(A^A)"'A^  (3.17) 

Pre-multiplying  a  vector  x  by  P.^  is  equivalent  to  projecting  the  vector  x  onto  the 
subspace  spanned  by  A.  Thus  the  result  x  =  P^x  lies  entirely  in  the  subspace  of  A 
and  the  difference  (or  error)  vector  (x  —  x)  lies  in  the  subspace  orthogonal  to  A,  see 
also  Remark  (ii)  above. 

3.4  Least  Squares  Solutions 

In  this  section  we  will  be  looking  at  the  deterministic  linear  equation  problem  of 


X  —  As  -I-  t; 


(3.18) 
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Given  that  x  and  A  is  known,  the  problem  is  to  estimate  s  such  that  the  squared 
error  Tr{{x  —  As)(x  —  As)^}  is  minimized  for  some  estimate  s. 

When  A  is  a  square  matrix  and  non-singular,  the  least  squares  solution  may  be 
found  by  inverting  A.  Thus, 

s  =  A"‘x.  (3.19) 

However,  when  A  is  not  a  square  matrix  but  is  rectangular,  then  the  least  squares 
solution  is  found  by  using  the  pseudo- inverse  of  A.  Thus, 

s  =  A*x.  (3.20) 

where  A*  is  the  pseudo-inverse  of  A.  The  pseudo- inverse  of  A  fulfills  the  following 
conditions,  see  [15,  51]. 

(I)  A* A  =  I 

(II)  AA*  =  Pa,  the  orthogonal  projector  onto  the  subspace  of  A. 

(III)  AA#A  =  A 

(IV)  A#AA*  =  A* 

When  A  has  rank  n  and  is  of  dimension  m-by-n,  (ie..  it  has  full  column  rank),  we 
can  write  the  pseudo- inverse  of  A  as 

A*  =  (AU)-‘A^  (3.21) 

Thus  the  least  squares  solution  for  s  as  found  from  Eq.  (3.20)  is  s  =  (A^A)"’ A^x. 

An  equivalent  formulation  of  the  least  squares  problem  would  be  to  find  the  closest 
estimate  x  which  lie  within  the  subspace  spanned  by  A  to  the  actual  observation  x. 
Then  the  solution  for  x  would  be  the  projection  of  the  actual  observation  x  onto  the 
subspace  of  A  ,  i.e., 


X  =  Pax 


24 


=  A(A’A)-‘A’x  (3.22) 

The  source  signal  estimate  corresponding  to  this  estimated  observation  is  observed 
to  be  the  same  as  in  Eq.(3.20). 


3.5  Matrix  Perturbation  Theory  Results 


In  this  section  we  will  state  first-order  perturbation  results  for  the  eigenvectors  and 
eigenvalues  and  pseudo-inverse  of  a  matrix  perturbed  from  a  nominal  matrix. 

3.5.1  Perturbation  of  Eigenvalues  and  Eigenvectors  of  a 
Hermitian  matrix 


Given  a  Hermitian  matrix  R  =  R  E  where  R  has  eigenvalues  {A*},  I  =  1,2,  ..,m 
and  eigenvectors  {ujt},/:  =  1,2,  and  the  perturbation  matrix  E  is  Hermitian 
and  small  enough,  see  [15],  then  the  eigenvalues  and  eigenvectors  of  the  perturbed 
matrix  can  be  written  to  first  order  as 


Afc  ~  Afc  -I-  u^Eut 

u|Euit 

Ufc  ^k  +  Y.-y — r 


(3.23) 

(3.24) 


In  practice,  at  least  in  statistical  signal  processing  applications,  the  covariance 
matrix  R  is  often  unknown  and  hence  have  to  be  estimated  from  finite  data.  Thus, 
assuming  that  the  data  is  wide-sense  stationary  and  that  there  are  N  available  samples 
or  snapshots  of  the  observation  array,  the  estimated  covariance  miatrix  is 


R=  ;^I^x(i)x(t)^ 


(3.25) 


ial 


When  the  observation  vector  x(i)  is  assumed  to  have  a  zero-mean  gaussian  distri¬ 
bution,  asymptotically  (for  large  N) 


E{K}  =  \i  +  0{N-^) 


(3.26) 
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Cov{X,X}  =  ^6.,  +  0(A^-')  (3.27) 

where  S,k  denotes  the  dirac  function  which  is  unity  for  i  —  k  and  zero  otherwise. 

E{ui}  =  u. +  0(iV-‘)  (3.28) 

Cot;{u.-,u*}  =  6,*:^f;_^i_u,uUO(iV-‘)  (3.29) 

See  [23]. 

3.5.2  Perturbations  in  Pseudo-Inverse 

Given  a  matrix  A  of  full  column  rank  n,  we  have  seen  in  the  t,revious  section  that 
its  pseudo-inverse  may  be  written  as  A*  =  (A^A)~‘A^  Suppose  that  instead  of  A. 
one  is  given  a  perturbed  matrix  A  =  A  -i-  AA.  Then  a  first-order  approximation  for 
the  pseudo-inverse  of  A  may  be  written,  see  [19],  as 

A’^  =  A*  +  (A^A)-‘ AA^I  -  AA’^)  -  A*AAA^  (3.30) 

3.6  Complex  Differentiation 

.■\rray  signal  processing  involves  the  computation  of  complex  quantitie.  Since  there 
will  be  a  need  to  take  the  derivatives  or  gradient  of  certain  function  with  respect  to 
complex  quantities,  this  section  presents  a  summary  of  the  operation,  see  [18]. 

Let  the  complex  quantity  z  be  written  as  2  =  +  jzi,  where  Zr,  Zi  denotes  the 

read  and  imaginaury  parts  of  z  respectively.  Then  the  derivative  operators  with  respect 
to  z  and  z'  are 


_  d  .  d 

(3.31) 

dz 

dZr  ^  ^  dzi 

d 

d  .  d 

(3.32) 

dz" 

dzy  ^  dzi 

This  operation  is  extended  to  the  vector  gradient  operation  in  the  following  manner. 
Suppose  that  w  is  a  complex  column  vector,  with  w  =  [wi  wj  •  •  ■  w^)*.  Then  the 
gradient  operator  is  defined  to  be  ^  =  [jJ;- 

Thus 


5(ctw) 

3w 

5(wtc) 

&w 

5(wtRw) 

dw 


0 

2c 

2Rw 


(3.33) 

(3.34) 

(3.35) 
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ASPECTS  OF  ARRAY  SIGNAL 
PROCESSING 

In  this  chapter  we  will  cover  ^  everal  conventional  nmethods  for  array  processing  such 
as  beamforming,  adaptive  beamforming,  eigenstructure  based  methods  such  as  MU¬ 
SIC  and  ESPRIT  and  maximum  likelihood  methods  (assuming  gaussian  distribution 
for  the  observation  vector).  Also,  the  estimation  of  the  number  of  sources  using  in¬ 
formation  theoretic  criteria  such  as  the  Akaike  Information  Criterion  (AIC)  and  the 
Minimum  Description  Length  (MDL)  will  be  discussed.  The  conventional  methods 
are  first  discussed  in  the  narrowband  model  case,  and  where  the  signal  sources  are  not 
roherent  with  one  another  (i.e.,  the  source  signal  covariance  matrix  is  non-singular). 
Extensions  of  these  methods  to  the  case  where  the  signal  sources  are  coherent  through 
the  use  of  spatial  averaging  using  sub-arrays  is  covered.  Finally,  the  extension  to  the 
wideband  case  is  also  discussed. 

4.1  Beamforming  Approaches 

Beamforming  is  a  spatial  filtering  approaches  to  array  processing.  Given  that  a  de¬ 
sired  source  signal  is  to  be  extracted  from  a  desired  DOA  angle,  the  outputs  of  the 
observation  array  are  delayed  and  summed  in  such  a  manner  as  to  direct  a  narrow 
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beam  of  sensitivity  for  signals  coming  from  the  desired  DOA  and  suppressing  others. 
Since  time  delays  may  be  equivalently  represented  as  phase  shifts,  the  output  of  a 
beam-former  may  also  be  viewed  as  the  result  of  an  inner  product  of  the  complex 
valued  beamformer  weight  vector  and  the  array  observation  vector,  as  in 

y  =  w^x  (4.1) 

where  the  column  vector  w  denotes  the  beamformer  weight  vector.  Note  that  the 
weight  vector  w  must  be  a  function  of  the  desired  DOA  angle.  It  may  or  may  not 
be  dependent  on  the  statistics  of  the  observation  vector  x  and  this  is  the  manner 
by  which  the  different  types  of  beamforming  approaches  may  be  classified,  see  [83j. 
Beamforming  approaches  where  the  weight  vector  is  selected  independent  of  the  ob¬ 
servations  are  classified  as  data  independent,  while  beamforming  approaches  whose 
weight  vectors  are  dependent  on  the  statistics  of  the  observation  data  are  classified 
as  statistically  optimum. 

4.1.1  Conventional  Beamformer 

The  conventional  beamformer  approach,  also  known  as  the  delay  and  sum  beamformer 
is  one  of  the  earliest  approach  used  for  array  signal  processing.  This  approach  is 
classified  in  the  data  independent  category  of  beamformers  as  the  beamformer  weights 
cire  not  dependent  on  the  observed  data  or  its  statistics. 

Intuitively,  the  idea  is  to  time  delay  the  different  array  sensor  outputs  and  sum 
them  in  such  a  manner  that  the  source  signal  emitted  from  the  desired  DOA  would 
add  coherently  while  those  of  the  interfering  and  noise  signals  would  add  incoherently. 
Thus,  the  resulting  beamformer  output  would  increase  the  SNR  of  the  desired  source 
signal  relative  to  the  other  signals,  including  noise. 

This  approach  can  be  illustrated  easily  for  a  ULA  with  a  single  source  signal,  at 
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DOA  angle  Assume  that  the  channel  gains  a,*  are  all  unity.  Thus,  the  z-th  sensor 
output  is 

x.(A:)  =  £-■'^<*-‘>^•’“’>5,  (4.2) 

Thus  if  the  beamformer  is  steered  towards  a  DOA  angle  of  0,  its  weight  vector  w 
is  defined  such  that 

W  =  [1.0  3^ 

Thus  the  output  of  the  beamformer  is 

y  =  W^x  =  Si  ^4  4^ 

Looking  at  the  power  of  the  beamformer  output,  we  have 

(4.5) 

Looking  at  the  above  expression  for  the  beamformer  output  power,  it  is  clear 
that  it  has  a  ’sinc-like’  oehavior  with  respect  to  the  argument  (sm^i  —  sinO).  The 
maximum  power  is  found  when  the  beamformer  is  steered  towards  the  source  signal 
at  the  DOA  0  =  0i,  and  that  it  has  a  main  lobe  width  of  The  conventional 
beamformer  would  therefore  not  be  able  to  resolve  a  case  there  are  two  source  signals 
with  DOA  angles  such  that  their  separation  is  less  than  half  the  main  lobe  width. 
Fig.  4.1  shows  an  example  of  a  sensitivity  plot  for  a  ULA  with  ten  sensors  spaced 
half  a  wavelength  apart  with  the  beam  steered  towards  30  degrees.  Note  that  any 
two  signals  with  DOA  angles  that  fall  within  the  main  lobe  would  not  be  resolvable. 
The  Rayleigh  resolution  limit  for  the  ULA  array  is  ^  (in  radians).  Note  that  this 
inability  of  the  conventional  beamformer  to  resolve  two  signals  with  close  DOA  angles 
is  related  to  the  array  aperture  size  and  not  due  to  the  amount  of  observation  data 
available. 


\  stn{^(sm5i  -  smd)} 
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Figure  4.1;  Conventional  Beamformer  Example 

When  a  general  array  geometry  is  assumed  instead  of  the  ULA  in  the  discussion 
above,  the  beamformer  weight  vector  is  defined  accordingly  as 

w  =  [1,0  (4.6) 

where  the  delays  r,  are  defined  as  in  Eq.(2.12). 

When  N  samples  or  snapshots  of  the  array  observations  are  available,  the  smoothed 
beamformer  power  estimate  is: 

m  = 

fesX 

=  wtRw  (4.7) 

where  R  =  jj  ELi  xx^{k)- 

The  conventional  beamformer  approach  is  thus  summarized  as  follows: 


1.  Estimate  R=  Efxx^}. 

2.  P{0)  =  w^Rw.  where  w(0)  =  [1  exp(-ju;oTi  ),••■.  exp{-ju)oTm-\  )]^- 

OOA  angle  estimates  are  found  from  peaks  of  P{6). 

Remarks: 

(i)  As  can  be  seen  from  Eq.(4.4),  the  conventional  beamformer  may  be  imple¬ 
mented  as  a  EFT  in  the  ULA  case.  Thus  it  is  an  attractive  approach  due  to  its 
simplicity  of  implementation. 

(ii)  Conventional  beamformers  suffer  from  the  problem  of  resolution,  i.e.,  source 
signals  that  are  situated  close  to  one  another  in  DOA  angles  cannot  be  resolved. 

(iii)  In  general,  conventional  beamformer  performance  is  degraded  by  the  presence 
of  multiple  sources.  It  performs  well  however  for  a  single  source  or  for  widely  separated 
sources. 

(iv)  ’Leakages’  or  spurious  signal  power  may  arise  due  to  the  relative*^  high  side- 
lobe  levels  of  the  conventional  beamformer.  The  sidelobe  levels  may  be  reduced  by 
using  ’window’  functions  (also  known  as  taper  weights)  in  the  beamformer  weight 
vector.  This  is  achieved  at  the  expense  of  reduced  resolution  capability  however  due 
to  increased  main  lobe  width. 

4.1.2  Linearly  Constrained  Minimum  Variance  Beamformer 

The  linearly  Constrained  Minimum  Variance  (LCMV)  beamforming  approach  is  based 
on  minimizing  the  beamformer  output  power  subject  to  certain  linear  constraints, 
[83].  Thus  interfering  signals  disjoint  from  the  source  signal  of  interest  would  have 
their  power  minimized.  These  constraints  may  consists  of  point  constraints  such  as 
constraining  the  signals  from  certain  DOA  angles  be  passed  through  the  beamformer 
with  specified  gains  and/or  phases,  or,  that  certain  inferference  signals  from  certain 
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DOA  angles  be  suppressed.  Derivative  constraints  may  also  be  used  to  control  the 
beamformer  response  over  a  specified  range  of  DOA  angles.  The  LCMV  beamformer 
problem  is  stated  as  follows. 

nnnw^Rw  such  that  C^w  =  f  (4.8) 

C  is  the  m-by-p  constraint  matrix  and  f  is  the  p-by-1  specified  response  vector 
(there  are  p  linear  contraints  on  the  problem).  Using  Lagrange  multipliers  for  the 
constraints,  we  write  the  minimization  problem  as 

4  =  w^Rw -f  (f  -  CW)U  (4.9) 

where  A  =  [Ai  Aj  •  •  •  Ap]‘  is  vector  of  Lagrange  multipliers. 

Taking  the  gradient  of  ^  with  respect  to  the  beamformer  weight  vector  and  setting 
it  to  zero  (a  necessary  condition  for  the  minimum),  we  get 

w  =  R-‘CA  (4.10) 

Using  the  constraint  relation  CW  =  f  and  assuming  that  the  constraint  matrix  C 
has  full  column  rank  and  the  covariance  matrix  R  is  invertible,  we  get 

Beamformer  Weight  Vector;  w  =  R“^C(C^R'’^C)~'f  (4.11) 

Beamformer  Power:  w^Rw  =  f^(C^R“^C)“*f  (4.12) 

It  is  clear  that  the  beamformer  weight  vector  w  is  dependent  on  the  covariance 
of  the  observation  vector  R  and  is  therefore  classified  in  the  statistically  optimum 
category  of  beamformers. 

The  LCMV  beamformer  weight  vector  w  may  be  decomposed  into  a  sum  of  two 
orthogonal  vectors  wc  and  W|i  which  are  their  projections  onto  the  subspace  spanned 
by  C  and  its  orthogonal  complement,  see  {83j.  From  Eq.(3.17)  in  the  previous  chapter. 
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we  have  that  the  orthogonal  projector  onto  the  subspace  spanned  by  the  constraint 
matrix  C  =  C(C^C)~'C^  and  thus 


Wc  =  PcW 

=  C{C^C)~‘f  (4.13) 

where  Wc  is  independent  of  the  observation  datal  Thus  the  only  component  of  w 
which  is  dependent  on  the  observation  data  statistics  is  =  P^w.  Since  is 
orthogonal  to  C,  we  may  write  it  as 

=  Db  (4.14) 

where  D  is  a  m-by-(m  -  p)  matrix  orthogonal  to  C  (it  can  be  found  by  applying  an 
SVD  decomposition  to  C)  and  b  is  a  (m  —  p)  column  vector.  Note  that  constrained 
minimization  of  Eq.(4.8)  is  reduced  to  the  un-constrained  minimization  as  below 

min(wc  +  Db)^R(w(7  +  Db)  (4.15) 

b 

since  C^w  =  C\wc  +  C^Db  =  f.  Solving  for  b,  we  get 

b  =  -(D^RD)-'D^Rwc  (4.16) 

Thus  the  LCMV  beamformer  may  be  implemented  with  a  combination  of  a  data 
independent  component  and  a  data  dependent  component  as  in  Fig.  4.2. 

This  equivalent  formulation  of  the  LCMV  beamformer  is  important  in  the  adaptive 
versions  of  the  LCMV  beamformer  as  will  be  shown  next. 

4.1.3  Adaptive  Beamforming 

Statistically  optimum  beamformers  as  in  the  LCMV  beamformer  above  require  an 
estimate  of  the  observation  covariance.  Since  the  statistics  of  the  data  may  vary  with 
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Dau  indqjendent  Part 


Data  Dq>endent  Part 


LCMV  Beamformer 

Figure  4.2:  Block  Diagram  of  LCMV  Beamformer 

time,  adaptive  beamformers  have  been  proposed  that  will  track  such  changes.  The 
problem  may  be  stated  as  follows,  see  (83): 

min4(w)  =  min£{l  z  -  w^u  I*} 

=  niin  +  w^IUw}  (4.17) 

Here  z  is  the  desired  respones,  u  is  an  input  data  vector  and  w  is  the  beamformer 
weight  vector.  is  the  power  of  the  desired  output,  =  E{uz*}  and  R,,  =  E{uu^}. 
The  solution  for  Eq.(4.17)  can  be  shown  to  be 

w  =  (4.18) 

Adaptive  methods  for  the  estimation  of  w  may  be  implemented  using  a  whole  host 
of  techniques  from  adaptive  filter  theory,  see  (11,  18).  Note  for  example  that  the 
Recursive  Least  Squares  method  or  the  Least  Mean  Squares  method  may  be  used, 
[18,  11). 

Using  the  Least  Mean  Squares  (LMS)  algorithm  to  allow  for  the  adaptive  updating 
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of  the  beamformer  weight  vector  estimate,  we  have 

Wjt+i  =  Wit  +  Qu{k)e'{k)  (4.19) 

where  e{k)  =  z{k)  —  w^(fc)u  and  a  is  the  step  size  or  gain  constant. 

Example  (Adaptive  LCMV  Beamformer):  The  LCMV  beamformer  consists 
of  a  data  independent  component  wc  which  may  be  computed  beforehand  and  a 
data  dependent  component  w^which  is  dependent  on  the  second  order  statistics  of 
the  data,  see  Eqs.(4.13),(4.16).  Thus  the  adaptive  version  of  the  LCMV  beamformer 
may  be  accomplished  by  defining  the  quantities 

Solution  Vector;  w  =  b 
Desired  Response:  z  =  w^x 
Input  Vector:  u  =  D^x 

where  the  quantities  b,  wc  and  D  are  as  defined  in  the  previous  section  on  LCMV 
beamformers. 

4.2  Eigenstructure  Based  Approaches 

In  the  past  two  decades  there  have  been  increasing  interest  in  eigenstructure  based 
approaches  for  array  signal  processing.  Unlike  the  conventional  beamforming  ap¬ 
proach  or  the  LCMV  beamforming  approach,  the  eigenstructure  based  approaches 
are  capable  of  (theoretically)  resolving  very  closely  separated  source  signals.  This  has 
led  to  these  methods  being  called  high  resolution  or  super  resolution  methods,  since 
they  are  capable  of  resolving  source  signals  with  DOA  separation  of  less  than  the 
Rayleigh  resolution  limit  as  found  in  a  conventional  beamformer. 

The  two  eigenstructure  based  methods  that  will  be  covered  in  this  section  are 
called  Multiple  Signal  Classification  (MUSIC)  and  Estimation  of  Signals  Parameters 


via  Rotational  Invariance  Techniques  (ESPRIT),  see  (52,  50|.  Before  we  proceed  to 
a  discussion  of  these  two  methods,  we  shall  first  discuss  the  basic  subspace  relations 
upon  which  these  two  methods  rest. 

Consider  the  narrowband  observation  covariance  as  in  Eq.  (3.8)  and  its  EED 
decomposition 

R  =  £;{xx^}  =  AQA^  + 

=  U,  A,U:  +  (4.20) 

Here  U,  of  dimension  m-by-n,U„  of  dimension  m-by-(m  —  n)  are  the  signal  subspace 
iind  noise  subspace  eigenvector  matrix  respectively.  The  signal  subsp2ice  eigenvector 
matrix  U,  spans  the  column  vectors  of  the  steering  matrix  A  and  thus  there  exists  a 
non-singular  matrix  T  such  that 

U,  =  AT  (4.21) 

This  is  easily  shown  by  considering  the  following.  Suppose  that  there  is  no  noise,  and 
thus  (7*  =  0.  Then  we  have 

U.A^Uj  =  AQA^  (4.22) 

Then  by  post  multiplying  both  sides  of  the  equation  above  with  U,A7’  >  we  get  Eq. 
(4.21),  where  T  =  QA^UjAj'.  Note  that  T  is  square  and  of  full  rank  (it  is  non¬ 
singular).  It  is  also  straightforward  show  that 

U.U^A  =  A  (4.23) 

It  is  also  clear  that  the  noise  subspace  eigenvector  matrix  U„  is  orthogonal  to  the 
steering  matrix  A,  i.e.,  U„  X  A,  since 


RU„  =  AQA^Un  4- (T^Un 


This  implies  that 

A»U„  =  0  (4.25) 

These  relations  comprise  the  backbone  of  most  eigenstructure  based  algorithms 
such  as  MUSIC  or  ESPRIT.  The  signal  subspace  relation  of  Eq.  (4.21)  is  used  in  the 
development  of  ESPRIT  while  the  orthogonality  of  the  noise  subspace  to  the  steering 
matrix  as  in  Eq.(4.25)  is  used  in  MUSIC. 


4.2.1  Multiple  Signal  Classification  (MUSIC) 


MUSIC  was  proposed  by  Schmidt  in,  see  (52,  53],  and  independently  by  Reddi.  see 
[45].  It  exploits  the  orthogonality  of  the  noise  subspace  to  the  steering  matrix  by  defin¬ 
ing  a  cost  which  measures  the  “closeness”  to  orthogonality  of  an  estimated  steering 
matrix  relative  to  the  noise  subspace.  From  Eq.(4.25),  it  is  clear  that 

Un^a(d)  =  0,  0  for  1  =  1, 2,  •  •  • ,  n  (4.26) 


where  a{^)  is  a  steering  vector. 

The  MUSIC  cost  function  is  thus  defined  to  be 


P{0)  = 


I 

II  uU(9)  r 


(4.27) 


When  the  array  manifold  is  known,  i.e.,  the  collection  of  all  possible  steering  vectors 
si(6)  for  all  possible  DOA  angles  is  known,  the  n  DOA  estimates  may  be  found  from 
the  n  highest  peaks  of  the  cost  function  P{9)  as  the  DOA  angle  parameter  is  swept. 
Thus  MUSIC  replaces  the  multidimensional  search  as  implied  in  Eq.(4.25)  with  a  one 
dimensional  search  over  the  range  of  possible  DOA  angles. 


Remark:  MUSIC  requires  the  array  manifold  to  be  known  or  measured  (cal¬ 
ibrated)  beforehand.  Thus,  while  it  has  the  attractive  feature  that  it  is  an  one¬ 
dimensional  search  procedure,  it  does  require  a  calibration  table  (a  lookup  table)  of 
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the  array  manifold  values  which  may  prove  costly  in  terms  of  memory  size  (depending 
on  the  grid  size  of  measured  calibration  values). 

In  theory  the  orthogonality  relation  in  Eq.(4.26)  should  hold  regardless  of  how 
close  the  angle  separation  of  the  DOA  of  any  pairs  of  source  signals  are.  Thus  MUSIC 
is  capable  of  resolving  source  signals  with  DOA  separation  that  are  much  smaller  than 
that  of  the  conventional  beamformer.  Under  finite  though  large  sample  conditions, 
however,  it  has  been  shown  that  MUSIC  does  exhibit  a  resolution  threshold  behavior, 
[23). 

In  general,  a  one  dimensional  search  is  necessary  in  the  implementation  of  MUSIC 
to  the  array  signal  processing  problem.  In  the  case  of  ULA  with  identical  sensors 
however,  as  may  have  been  noticed  by  the  reader  by  now,  an  elegant  method  may  be 
used  instead  of  searching  for  the  DOA  angles.  In  this  case,  a  variant  of  MUSIC  called 
Root-MUSIC  can  be  applied.  Note  that  the  steering  vector  a(z)(0)  is  Vandermonde 
and  may  be  written  as  a(0)  =  [1.0  2*  •  •  •  where  z  s  Then 

Eq.(4.26)  may  be  written  as 

U].a(2)  =  0,  at  r  =  (4.28) 

Form  the  polynomial 

D(.-)  =11  Ul,a(r)  IP  (4.29) 

The  polynomial  /?(z)  =  0  for  2  =  i.e.,  it  has  zeros  on  the  unit  circle 

corresponding  to  the  true  DOA  angles.  Thus  the  Root-MUSIC  method,  see  [3],  forms 
the  polynomial  D(z)  from  the  estimate  of  the  noise  subspace  eigenvector  matrix  and 
solves  for  roots  lying  on  the  unit  circle. 

In  the  past  few  years,  the  performance  analyses  of  MUSIC  and  its  variants  have 
been  a  subject  of  intense  interest.  Of  particular  note  is  the  work  found  in  [23],  which 
is  one  of  the  first  papers  to  analyze  MUSIC.  See  also  [62,  63,  43],  What  is  clear  from 
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these  analyses  is  that  MUSIC  provides  DOA  estimates  that  are  asymptotically  un¬ 
biased,  gaussian  distributed  and  in  general  is  not  equivalent  to  Maximum-Likelihood 
estimates,  see  next  section  and  also  {64,  62).  In  the  case  of  uncorrelated  source 
signals  however,  MUSIC  and  ML  methods  are  asymptotically  equal  in  MSE.  The 
Root-MUSIC  and  MUSIC  methods  are  asymptotically  equivalent  as  has  been  shown 
in  [43]. 


4.2.2  Estimation  of  Signal  Parameters  via  Rotational  In 
variance  Techniques  (ESPRIT) 


In  contrast  to  MUSIC  which  utilizes  the  noise  subspace  relation  of  Eq.(4.25),  ESPRIT 
exploits  the  signal  subspace  relation  of  Eq.(4.21),  namely  the  fact  that  the  steering 
matrix  A  is  spanned  by  the  signal  subspace  eigenvector  matrix.  The  rotational  in¬ 
variance  component  comes  into  play  here  because  ESPRIT  tissumes  that  A  has  a 
certain  special  invariance  structure,  see  {50].  Specifically,  A  must  have  the  structure 
that  it  may  be  partitioned  in  the  following  manner. 


A  = 


B 

m 


(4.30) 


where  B  is  a  steering  sub-matrix  and  $  is  a  diagonal  matrix  such  that 


(4.31) 


L  O  J 

Here  ESPRIT  assumes  an  array  structure  composed  of  two  sub-arrays  we  shall  call 
X  and  Y  sub-arrays,  where  each  sensor  in  X  is  matched  (identical)  to  its  pair  in  Y 
and  spaced  at  a  distance  d  from  each  other,  see  Fig.  4.3  for  example.  Note  that  only 
the  sensor  pairs  (also  known  as  doublets)  need  to  be  matched  and  different  sensor 
characteristics  for  different  pairs  are  allowed.  Using  simple  geometry  to  calculate 
the  time  delay  difference  between  the  sensors  of  a  sensor  pair,  we  can  see  that  the 
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Figure  4.3:  Esprit  Subarrays  Example 

difference  in  time  delay  for  a  signal  coming  from  DOA  angle  6k  at  the  t-th  sensor  in  X 
and  its  pair  in  Y  is  At,*  =  Thus  if  the  observation  column  vector  is  arranged 

by  concatenating  the  observations  from  the  X  sub-array  on  top  of  the  observations 
from  the  Y  sub-array,  then  the  steering  matrix  A  which  results  will  fulfill  the  required 
rotational  invariance  as  in  Eq.(4.30).  From  the  fact  that  the  steering  matrix  A  is 
spanned  by  the  signal  subspace  eigenvector  matrix,  see  Eq.  (4.21),  there  exists  a 
non-singular  matrix  T  such  that 

AT  = 


where  ExiEv  are  matrices  partitioned  from  U,  with  dimensions  corresponding  to 
the  size  of  the  X  and  Y  sub-matrices.  Furthermore,  the  matrices  Ex,Ev  both  span 


U, 


B 

Ex 

Ey 


(4.32) 
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the  same  subspace  due  to  the  structure  of  the  steering  matrix  A.  Thus,  there  exists 
a  matrix  such  that 

Ex9  =  Ey  (4.33) 

Notice  however  that  due  from  Eq.(4.32) 

BT«’  =  B$T  (4.34) 

Since  T  is  non-singular  and  B  has  full  column  rank  this  implies  that  the  eigenvalues 
of  ^  are  exactly  the  diagonal  terns  in  $!  Thus  the  ESPRIT  algorithm  comprises 
primarily  the  estimation  of  the  signal  subspace  eigenvector,  its  partitioning  into  sub- 
matrices  E.v,Ev  solving  for  ’i'  from  Eq.(4.33)  and  finally  the  extraction  of  the  eigen¬ 
values  of  9.  The  DO  A  angles  may  thus  be  extracted  from  the  phase  terms  in  the 
eigenvalues. 

Remarks: 

(i)  may  be  solved  from  Eq.(4.33)  by  using  either  Total  Least  Squares  (see 
[50,  15])  in  which  case  the  algorithm  is  called  TLS-ESPRIT  or  by  ordinary  least 
squares  in  which  Ccise  it  is  called  LS-ESPRIT. 

(ii)  ESPRIT  circumvents  the  need  for  a  search  procedure  and  for  knowledge  of /or 
calibration  of  the  array  manifold.  Thus  it  is  an  attractive  alternative  to  MUSIC.  It 
requires  however  the  special  rotational  invariance  structure  of  the  array  which  entails 
matched  sensor  pairs  (matched  in  terms  of  characteristics  and  alignments). 

Performance  analysis  of  ESPRIT  h.i.ve  been  conducted  by  Stoica  et.,  al,  [63],  Rao 
et.,  al,  [44],  Ottersten  et.,  al,  [35].  One  interesting  result  of  Stoica/Nehorai  is  that 
the  performance  of  MUSIC  is  aymptotically  better  than  that  of  ESPRIT.  The  DOA 
estimates  have  been  shown  to  be  zero  mean,  gaussian  distibuted,  see  for  example 
[35].  It  has  also  been  shown  that  asymptotically,  TLS-ESPRIT  and  LS-ESPRIT  are 
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equivalent,  although  simulations  have  suggests  that  TLS-ESPRIT  is  superior  in  low 
SNR  and  small  snapshot  size  cases,  see  [50,  44]. 

4.3  Maximum  Likelihood  (ML)  Methods 

There  are  two  formulations  for  maximum  likelihood  methods  for  the  narrowband  array 
signal  processing  problem,  the  so  called  stochastic  maximum  likelihood  formulation 
and  the  deterministic  maximum  likelihood  formulation.  The  deterministic  tag  here 
refers  to  the  assumption  that  the  source  signals  themselves  are  deterministic,  the 
observations  however  are  still  random  since  the  additive  noise  is  AWGN.  Maximum 
Likelihood  methods  are  of  interest  in  all  facets  of  estimation  theory  since  under  certain 
conditions,  ML  estimates  have  mean  square  errors  that  are  asymptotically  equal  to 
the  Cramer- Rao  lower  bounds,  see  [29], 

4.3.1  Stochastic  Maximum  Likelihood 

In  the  stochastic  ML  formulation  the  source  signals  have  zero  mean,  temporally  white 
gaussian  distributions  such  that 

=  QV  (4.35) 

£:{5(()s'(0)  =  0  (4.36) 

Writing  the  probability  density  function  for  N  snapshots  of  the  observation  array,  we 
get 

/©(xi,---,x^)  = -  (4.37) 

(27r)/"/ydct(R)  I  2,tl  '  J 

where  ©  is  a  column  vector  of  the  parameters  to  be  estimated  which  includes  the  n 
DOA  angles,  the  entries  of  the  source  covariance  matrix  Q  and  the  noise  variance  <t^. 
The  ma.’:imum  likelihood  method  entails  finding  the  estimate  of  the  parameter  vector 
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0  such  that  the  probability  density  function  in  Eqd-l.S?)  is  maximized.  Taking  the 
negative  log  of  the  pdf  (since  the  log  function  is  monotonic)  this  is  equivalent  to 


,v 


0  =  min 

0 


log{det{R))  +  xIR  ‘xjt 


ibsl 


(4.38) 


Allowing  for  the  fact  that  the  source  covariance  matrix  is  Hermitian,  the  total  num¬ 
ber  of  parameters  that  must  be  estimated  is  n  DOA  angles+m^  real  and  imaginary 
parts  of  source  covariance  matrix -fl  noise  variance.  This  minimization  is  thus  often 
computationaly  expensive  as  it  essentially  means  that  a  multi-dimensional  (specifi¬ 
cally  (n  +  -f  1))  search  must  be  conducted.  The  non-linearity  of  the  minimization 
criterion  of  Eq.(4.38)  further  compounds  the  problem  due  to  the  likely  existence  of 
local  minima  solutions. 


For  methods  proposed  to  solve  for  the  estimates  of  the  stochastic  ML  estimates, 
see  [34]  and  the  references  therein.  For  performance  analysis,  see  [34]  -  the  stochastic 
ML  method  is  asymptotically  efficient  if  the  global  minimum  of  Eq.(4.38)  is  found. 
See  also  [64,  62]. 


4.3.2  Deterministic  Maximum  Likelihood 


In  the  deterministic  ML  formulation,  the  source  signals  are  assumed  to  be  determin¬ 
istic  signals,  i.e.,  they  are  not  random  processes.  The  observations  however  would 
still  be  gaussian  distributed  due  to  the  additive  white  gaussian  noise  component.  As¬ 
suming  that  there  are  N  snapshots  of  the  observation  array  vector,  we  can  write  the 
probability  density  function  of  the  observations  as 

/e(xt.---,x/v)  =  i^exp|-^^^(x,  -  Asfc)^(xfc  -  Asfe)|  (4.39) 

where  Sk,k  =  1,  •  ■  • ,  jV  are  the  deterministic  source  signal  vectors  and  A  is  the  steering 
matrix.  The  parameters  within  the  parameter  vector  0  that  must  be  estimated  are 
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the  n  DOA  angles,  the  N  signal  source  vectors  and  the  noise  variance,  a‘.  The 
of  Eq.(4.39)  may  also  be  viewed  as  a  conditional  pdf  (conditioned  on  N  realization 
of  the  source  signal  vector).  Taking  the  negative  log  of  Eq.(4.39),  we  get 

=  min  |miVlog<7^  +  -^  J^(xjt-ASit)^(Xfc- As*)]  (4.40) 

Taking  the  derivative  with  respect  to  the  noise  variance  cr^  and  setting  to  zero,  we 
can  solve  for  its  ML  estimate  as 

<T^  =  X^(xfc  -  Asfc)^(xit  -  Asfc)  (4.41) 

fc=i 

Next,  taking  the  derivative  of  the  log  likelihood  function  with  respect  to  the  deter¬ 
ministic  source  signal  vector  Sk  and  setting  to  zero,  find  that  their  ML  estimates 
are 

Sfc  =  (AU)-*A^x  (4.42) 

Plugging  these  ML  estimates  for  the  noise  variance  and  source  signal  vectors  into  the 
minmization  of  Eq.(4.40),  we  get  the  compressed  deterministic  ML  problem 

(«,)  =  rninTrlPiR}  (4.43) 

where  P;^  =  I  -  P>i  and  R-  =  ^  llk=i  Xtx[.. 

Remark:  The  deterministic  ML  problem  is  a  nonlinear  minimization  problem 
and  thus  computationally  expensive.  Various  methods  have  been  proposed  for  this 
problem,  see  [4,  84].  Statistical  performance  analysis  have  been  conducted,  see  [64, 
62,  34]  etc. 

4.4  Cramer  Rao  (CR)  Bounds 

Given  any  unbiased  estimator  of  the  parameter  vector  0,  the  error  covariance  ma¬ 
trix  of  its  estimated  parameter  vector  is  lower  bounded  by  the  Cramer  Rao  bound. 


Specifically,  given  any  estimate  ©  such  that 


E{®}  =  © 

andC  =  £’{(©-©)(©-©)'} 

then  the  error  covariance  matrix  C  is  lower  bounded  such  that 

C>J-'  (4.44) 

where  J  is  the  so-called  Fisher  Information  matrix  and  the  inequality  is  taken  to  mean 
that  C  — J“*  >  0,  i.e.,  it  is  non-negative  definite.  Let  the  pdf  of  the  observation  vector 
be  denoted  by  /©(x).  Then  the  Fisher  Information  matrix  is  defined  as 

The  CR  bound  thus  provides  a  benchmark  by  which  the  ’goodness’  of  an  unbiased 
estimator  may  be  measured.  This  is  why  maximum  likelihood  estimators  are  often 
of  interest  since  under  certain  regularity  conditions,  the  error  covariance  matrix  of 
a  maximuM  likelihood  estimator  will  approach  the  CR  bound  asymptotically,  i.e.,  it 
is  asymptotically  efficient.  An  unbiased  estimator  whose  error  covariance  matrix  is 
asymptotically  equal  to  the  CR  bound  is  an  asymptotically  efficient  estimator. 

In  general  the  evaluation  of  the  CR  bound  for  an  arbitrary  pdf  of  the  observation 
data  may  prove  to  be  a  formidable  task.  Fortunately,  for  the  case  where  the  pdf  is 
gaussian,  the  CR  bound  attains  relatively  simple  forms. 

We  can  consider  two  formulations  for  the  aissumption  of  gaussian  distributed  ob¬ 
servation  data  vectors,  namely  the  stochastic  where  the  source  signals  are  also  zero 
mean,  gaussian  distributed  (temporally  white)  as  discussed  in  the  previous  section 
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on  the  Stochastic  ML  estimator  and  the  deterministic  where  the  source  signals  are 
considered  to  be  deterministic,  as  discussed  in  the  section  on  the  Deterministic  ML 
estimator. 

4.4.1  Stochastic  CR  Bound 


Using  the  formulation  as  found  in  the  Stochastic  ML.  the  Fisher  Information  matrix 
J  can  be  shown  to  be  (see  [2,  51])  comprised  of  the  following  element  entries 


J.k 


a©,  aej 


(4.47) 


Here  N  denotes  the  number  of  snapshots  available  and  R  is  the  observation  covariance 
matrix. 


Remark:  For  explicit  evaluations  of  the  stochastic  ML  see  [64,  34)  etc. 


4.4.2  Deterministic  CR  Bound 


In  the  deterministic  case  with  the  formulation  as  in  the  Deterministic  ML  discussion 
previously,  the  CR  bound  matrix  have  been  evaluated  by  Stoica  and  Nehorai  in  [62] 
where 

_2  r  Af  "I 

J-'  =  T  [X^(l-)D^PiDX()l-)]  I  (4.48) 

where  X(fc)  =  Diag{x{k)},  i.e.,  a  diagonal  matrix  with  its  diagonal  entries  taken 
from  the  elements  of  the  observation  vector  x(k)  and  D  defined  such  that  its  i-th 
column  vector  is  d,  =  evaluated  at  0  =  6i. 

Remark:  Stoica  and  Nehorai  derived  the  above  expression  and  discussed  several 
important  special  cases  of  interest,  [62]. 
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4.5  Information  Theoretic  Model  Order  Identifi¬ 
cation 

Throughout  this  thesis  it  shall  be  assumed  that  the  number  of  source  signals  is  known, 
see  assumption  A4.  In  practice  however,  the  number  of  source  signals  would  also 
have  to  be  estimated  from  the  observed  data.  Two  model  order  estimation  criteria 
which  have  gained  a  lot  of  attention  in  the  literature  are  the  Akaike’s  Information 
Criterion  (AIC),  see  Akaike,  and  the  Minimum  Description  Length  (MDL)  criterion, 
see  Risscinnen,  Schwartz. 

Both  methods  attempt  to  find  the  model  order  that  best  fits  the  given  family  of 
pareuneterized  model  pdf  /(x  j  0).  AIC  is  proposed  as  a  measure  (see  .Akaike,  (1|) 
where  it  was  shown  that  minimizing  the  AIC  is  asymptoticiilly  equivalent  to  mini¬ 
mizing  the  Kullback-Leibler  function  which  is  a  measure  of  the  statistical  “distance” 
of  the  true  model  pdf  and  the  estimated  model  pdf.  MDL  (see  Rissanen,  Schwartz, 
[46,  54])  is  derived  by  treating  the  problem  as  that  of  determining  the  minimum  code 
length  that  will  describe  the  observation  data.  Thus  the  AIC  and  MDL  are  stated  as 

AlC(fc)  =  -2logf{x  1  -h  2k  (4.49) 

MDL(^')  =  -log fix  I  ©<*’)  -f-  hlogN  (4.50) 

where  ©  is  the  model  parameter  vector  and  ©**^  is  its  ML  estimate,  k  is  the  number 
of  free  adjusted  parameters  in  the  model  and  N  is  the  number  of  snapshots  of  the 
observation  data.  The  model  with  the  minimum  AIC  or  MDL  is  selected  as  the  one 
that  best  fits  the  observed  data. 

Wax  and  Kailath  [77]  had  applied  AIC  and  MDL  to  the  problem  of  determining  the 
number  of  source  signals  in  the  narrowband  array  signal  processing  problem.  Using 
the  statistiv-Ai  assumptions  in  the  Stochastic  ML  model,  these  criteria  are  shown  to 


AlC(n)  =  -2log 


nm  — n) 

l=rn-f  I  ^1 _ 

-i— V"*  .  A 

m— n  ^t=n+l  ’ 


|m— n)A^ 


4-  2n(2m  —  n) 


rm  ”) 


(m-n)N 


+  “n(2m  —  h)logN 


(4.52) 


MDL(ft)  =  -%  +ln{2m-h]logN  (4.52) 

\  m-n  ^«=n+l  } 

where  n  denotes  the  estimated  number  of  source  signals  and  A,  denotes  the  eigen 
values  of  the  estimated  covariance  matrix  of  the  array  observation  vector  R  arranged 
in  order  of  decreasing  magnitude. 

Thus  the  estimated  number  of  sources  are  picked  from  the  value  that  minimizes 


these  equations  depending  on  which  criterion  is  to  be  used. 


Remark:  Wax  and  Kailath,  [77]  showed  that  for  the  AIC  and  MDL  crueria  as 
in  Eqs.(4.51),  (4.52),  that  AIC  is  inconsistent,  ie.  asymptotically  it  tends  to  over¬ 
estimate  the  number  of  sources.  MDL  however  is  asymptotically  consistent.  Exten¬ 
sions  of  these  information  theoretic  model  order  estimation  methods  for  wideband, 
coherent  sources  may  be  found  in  [20,  78,  75]. 


4.6  Spatial  Averaging  for  Coherent  Sources 


In  the  development  of  the  eigenstructure  based  approaches  thus  far,  it  has  been 
assumed  that  no  two  source  signals  are  fully  coherent  with  each  other,  i.e.,  the  source 
signal  covariance  matrix  is  positive  definite  (assumption  A2).  In  the  situation  where 
two  or  more  source  signals  are  coherent,  the  application  of  MUSIC  would  only  resolve 
the  non-coherent  source  signals.  Similarly,  the  ESPRIT  algorithm  would  also  fail  as 
the  signal  subspace  eigenvector  matrix  would  have  a  rank  less  than  the  rank  of  the 
true  steering  matrix. 

In  the  case  of  coherent  sources,  the  technique  of  spatial  averaging  has  been  pro¬ 
posed  in  [56].  Spatial  averaging  assumes  the  use  of  a  ULA  with  m  identical  sensors 
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divided  into  L  overlapping  subarrays  each  of  size  p  sensors.  The  multiple  covariance 
matrices  of  these  sub-arrays  are  then  summed  and  divided  by  L,  thus  getting  an  av¬ 
erage  covariance  matrix  upon  which  the  eigenstructure  based  method  is  applied.  In 
[56]  it  was  shown  that  this  pre-processing  of  the  data  would  lead  to  the  averaged 
source  covariance  matrix  being  of  full  rank,  provided  that  the  number  of  subarrays 
are  greater  than  or  equal  to  the  number  of  sources,  L  >  n.  Note  that  the  number 
of  sensor  within  a  subarray  p  must  also  be  greater  than  or  equal  to  the  number  of 
sources  n.  Thus  the  use  of  spatial  averaging  to  alleviate  source  coherence  problems 
lead  to  a  tradeoff  with  respect  to  the  available  array  aperture  (i.e.,  it  leads  to  a  smaller 
aperture  being  used,  less  than  the  available  one)  since  L  subarray  covariance  matrices 
need  to  be  averaged.  It  has  been  shown  also  that  spatial  averaging  requires  at  least 
50  percent  more  sensors  for  the  same  effective  array  aperture  as  before,  see  [4l|. 

4.7  Wideband  Extensions 

Wideband  extensions  for  array  signal  processing  fall  roughly  into  two  camps  or  types 
of  approaches,  the  so-called  incoherent  and  coherent  approaches.  Both  are  applied  in 
the  frequency  domain  and  may  be  used  to  extend  most  eigenstructure  based  methods 
to  the  wideband  case. 

4.7.1  Incoherent  Approach 

Intuitively,  the  incoherent  approach  for  wideband  array  processing  may  be  summa¬ 
rized  as  follows:  take  the  frequency  transform  of  the  observation  data  which  divides 
the  data  into  several  frequency  bins.  Apply  the  desired  narrowband  estimation  tech¬ 
nique  then  on  data  at  various  frequency  bins.  Finally,  use  the  estimates  found  at  the 
different  frequency  bins  to  arrive  at  some  sort  of  average  estimate.  Thus  individual 
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narrowband  processing  is  done  at  various  frequency  bins  with  a  rinal  averaging  of  the 
results. 

Consider  an  extunple  of  this  wideband  extension  approach  for  beamfoming  appli¬ 
cations.  Note  that  first  an  FFT  is  taken  of  the  observation  data  and  then  individual 
narrowband  beamformers  are  applied  to  the  frequency  transformed  data.  Finally  the 
various  beaunformer  outputs  are  added  and  inverse  FFT  is  performed  yielding  the 
final  output  estimate  y(t). 

In  the  case  of  MUSIC,  an  incoherent  approach  was  proposed  in  [78],  where  the 
cost  function  P{d)  is  defined  as 

pm  =  — i - ! -  (4.53) 

ELi  II  U.(u,)ta(u,,9)  r 

Here  U„(w/t)  denotes  the  noise  subspace  eigenvector  matrix  of  R(a;fc),  the  observation 
covariance  matrix  at  frequency  Uk,  see  Eq.  (3.10). 

See  (65)  and  [78]  for  extensions  of  this  type  for  MUSIC.  Extensions  for  ESPRIT 
are  found  in  [33]. 

4.7.2  Coherent  Approach 

The  coherent  approach  for  wideband  processing  is  based  on  the  idea  that  the  obser¬ 
vation  data  at  the  various  frequeny  bins  should  be  added  coherently  in  some  manner 
before  the  narrowband  eigenstructure  based  method  is  applied.  Intuitively  the  addi¬ 
tion  of  the  information  at  the  various  frequency  bins  vould  tend  to  improve  the  SNR 
of  the  data.  This  approach  was  proposed  in  Wang  and  Kaveh,  see  [75]  . 

Define  the  transformation  matrices  Tfu;*)  such  that: 


T{uk)A{uk)  =  A{u;o),k  =  1.  •  •  • ,  J 


(4.54) 
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where  wq  denotes  the  focusing  frequency  to  which  the  steering  matrices  of  other 
frequencies  are  focused  to.  Note  that  this  requires  knowledge  of  the  true  DOA  angles 
which  is  exactly  the  quantities  to  be  estimated!  This  dilemma  is  circumvented  by 
using  a  rough  estimate  of  the  DOA  angles  such  as  that  found  by  using  conventional 
beamforming.  Then  the  various  covariance  matrices  at  the  different  frequencies  are 
summed  in  the  following  manner. 

R,  =  f;  T(u;fc)RK)T(a;*)^  (4.55) 

/c=l 

Finally  the  desired  narrowband  method  is  applied  on  the  focused  covariance  ma¬ 
trix  Rx-  Extensions  to  MUSIC  for  this  approach  can  be  found  in  [20.  S.  27].  The 
extensions  for  ESPRIT  can  be  found  in  [21]. 
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ANALYSIS  OF  ESPRIT  UNDER 
RANDOM  MODEL  ERRORS 


5.1  Introduction 

£/igenstructure  baaed  algorithms  for  array  signal  processing  has  gained  considerable 
interest  in  recent  years.  Algorithms  such  as  MUSIC  [53]  and  ESPRIT  [50]  have  been 
proposed  generally  for  the  DOA  estimation  problem.  Much  work  on  performance 
suialysis  of  these  algorithms  has  been  directed  to  studying  the  effects  of  noise  and 
finite  number  of  snapshots,  with  no  errors  assumed  on  the  sensors  gains/phases, 
locations  etc.  ,  see,  e.g.,  [23], [35]. 

This  chapter  presents  analysis  results  on  the  MSE  of  ESPRIT  DOA  angle  esti¬ 
mates  with  errors  in  the  array  sensors.  Such  errors  include  random  errors  in  the  sensor 
gains  and  phases,  random  errors  in  sensor  locations  and  random  errors  in  sensor  pair 
alignments.  Both  MUSIC  and  ESPRIT  perform  well  when  the  sensors  in  the  array 
are  either  calibrated  (i.e.,  the  array  manifold  is  known)  as  in  MUSIC  or  the  sensors  in 
a  sensor  pair  are  matched  perfectly  as  in  ESPRIT.  In  practice,  however,  this  may  not 
always  be  possible  due  perhaps  to  external  environmental  effects  on  the  array,  dete¬ 
rioration  of  electronic  components,  measurement  errors,  etc.  As  a  consequence,  there 
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has  been  interest  in  the  performance  analysis  of  these  methods  under  model  errors. 
Analysis  of  MUSIC  under  various  model  errors  can  be  found  in  {82]-[68).  Analysis 
of  ESPRIT  where  the  steering  matrix  errors  are  modeled  as  uncorrelated,  zero-mean 
Gaussian  distributed  can  be  found  in  [30,  68].  In  [68,  69],  the  analysis  of  ESPRIT 
under  model  errors  are  performed  using  its  formulation  as  the  minimization  of  a  cer¬ 
tain  cost  function.  The  MSE  expressions  obtained  therein  are  very  complicated  due 
to  its  parameterization  in  terms  of  not  only  the  DOA  angles  but  also  other  ’nuisance’ 
parameters  [35],  such  as  the  real  and  imaginary  parts  of  the  steering  matrix,  the  mag¬ 
nitudes  of  the  elements  of  a  diagonal  matrix,  etc.  In  this  chapter,  a  direct  approach 
towards  the  analysis  of  ESPRIT  under  model  errors  is  taken  where  the  DOA  angles 
are  the  sole  parameters  considered.  This  approach  yields  simple  MSE  expressions 
that  are  related  explicitly  to  the  model  errors.  As  such,  it  provides  interesting  insight 
into  the  performance  sensitivity  of  ESPRIT  to  model  errors. 

5.2  Problem  Formulation 

Throughout  this  chapter,  the  superscripts  ’H’,  T’,  and  denote  the  Hermitian, 
transpose,  complex  conjugate  and  the  Moore-Penrose  pseudo-inverse  of  a  matrix, 
respectively.  Also,  n  denotes  the  number  of  narrowband  far-held  sources,  m  denotes 
the  number  of  sensors  within  a  subarray,  mo  denotes  the  total  number  of  sensors 
(i.e.  mo  <  2m)  and  {/x,,  i',}  and  {/I,,  t/,}  denote  the  sensor  locations  in  the  X  and  Y 
(possibly  overlapping)  subarrays  respectively. 

Here,  we  will  be  concerned  with  the  situation  where  the  covariance  matrix  is  as¬ 
sumed  known  (i.e.,  the  hnite  sample  effects  on  data  covariance  matrix  are  not  in¬ 
vestigated  here).  Under  nominal  conditions,  the  sensor  pairs  within  the  X  and  Y 
subarrays  are  matched  perfectly,  with  a  distance  spacing  d  between  them  [50].  By 
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appropriate  choice  of  coordinates,  one  can  rewrite  the  sensor  locations  within  the  X 
and  Y  subarrays  to  be  v,}  and  {/i,  +  respectively.  Then,  the  observation 

vectors  from  the  X  and  Y  subarrays  can  be  written  as: 

X  subarray  :  x  =  FAs  +  tjj.  and  Y  subarray  :  y  =  FA^s  +  rjy  (5.1) 

Here  F  denotes  the  diagonal  array  of  nominal  sensor  gains  and  phases,  namely, 

r  =  Diag{aie^*‘ , ..,  and  A  is  the  steering  matrix  defined  by  the  set  of  distinct 

DOA  angles  (the  DOA  angles  are  measured  with  reference  to  the  normal  of  the 
sensor  pairs), 

A  =  aj  ftn  j  .  (5.2) 

The  individual  steering  vector  at  DOA  angle  Ok  is  defined  as 

at  =  [exp{-j^{fiisin6k  +  I'lCosOk)}... 

27r 

•  •  •  exp{-j—{fim3in6k  +  iymCos6k)}f .  (5.3) 

Furthermore,  #  =  Diag{cxp(— j^sm^i),  ..,exp(— j^dsin0n)}»  ^  is  the  wavelength 
of  the  narrowband  signals,  s  is  the  source  vector  and  denote  the  additive  noise 
vectors  at  the  X  and  Y  subarrays,  respectively. 

It  is  assumed  that  the  matrix  A  has  full  column  rank,  the  narrowband  sources  are 
incoherent  and  the  noise  vectors  are  uncorrelated  zero  mean  processes  with  covariance 
see  assumptions  A1-A4. 

5.2.1  General  Formulation  of  Error  Models 

Under  small  array  model  errors  and  using  first  order  approximation,  a  general  for¬ 
mulation  of  the  array  signal  processing  problem  under  model  errors  can  be  derived, 
namely, 

X  subarray  :  x  =  (FA  -f  dXx)s-^T]x  and  Y  subarray  :  y  =  (FA  +  5Ay)#s+7/j,  (5.4) 


00 


The  5Ax,  dAy  terms  are  dependent  on  the  types  of  error  model  in  question.  These 
different  error  scenarios  are  derived  as  below. 

Errors  in  Sensor  Gains  and  Phases:  dA^  =  AFxA  and  dAy  =  AF^A  where 

AFx  =  DiagiiAax  +  joj  A<^,)e-'‘*' ,  •  •  • ,  (Aa„,  +  ja„A4>m)e^*”'}.  (5.5) 

AFy  =  Dzag{(Aat  -h  ;QiA<^i)e-'^'.  •  •  •  ,{AQ:m  -h  jamA^m)e^*”'}-  (5.6) 

Here  {(Aq,,  A^i)}  and  {(Aq,,  A^,)}  denote  errors  in  the  sensor  gains  and  phases  of 
the  X  and  Y  subarrays  respectively. 

proof:  Consider  the  f-th  perturbed  sensor  gains  and  phases:  (q,  + 

Using  first  order  approximation  with  ~  €■''”•(1  +  jA(/>i)  and  discarding  ail 

terms  higher  than  first  order,  the  error  terms  of  the  above  is  found. 

Using  the  same  approach  for  the  other  error  models,  we  get 
Errors  in  Sensor  Locations:  dA^  =  FAAx  and  ^Ay  =  FAAy  where 


2ir 

AAr  =  -J-j~(Dtag{A^i,- 

•  • ,  A^,n}A,  +  Diag{  Auir  •  •,Ai/m}Ac) . 

(5.7) 

2zr 

AAy  =  -;Y(Z)ia^{A/ii,- 

■  ’,A/i,n}As  +  Diag{  APi,-  ■  • .  Ai/m}Ac)  ■ 

(5.8) 

Here  {(A/i,,  Ai^,)}  and  {(A/i,,  Ai/i)}  denote  errors  in  the  sensor  locations  of  the  X  and 
Y  subarrays  respectively,  A,  =  [aisindi,...,ansm0„j  and  Ac  =  [ajcos^i, ...,anCOsd„]. 

Errors  in  Sensor  Pair  Alignments:  In  this  case,  the  sensor  pairs  within  the  X 
and  Y  subarrays  are  perturbed  by  small  alignment  errors  such  that  they  are  rotated 
by  small  angles  from  the  nominal  aligned  direction.  Hence,  dAx  =  FAAr  and  dAy 
=  FAAy  where 


A  Ax  =  -j-^Diag{^i,-  ■ 

•./3m}  Ac. 

(5.9) 

ird 

AAy  =  ~^j—Diag{f3u-- 

'  >  i^m  }  Ac* 

(5.10) 
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Here  {/?,}  denote  the  small  angle  rotations  from  the  normal  of  the  sensor  pairs  due 
to  misalignment  errors  and  Ac  is  as  defined  previously. 

Remark:  Note  that  the  random  sensor  error  assumptions  here  enter  directly  into 
the  sensor  parameters  such  as  gains/ phases,  locations  and  alignments  as  opposed  to 
the  uncorreiated  Gaussian  error  assumptions  imposed  on  the  steering  matrix  (i.e.,  on 
dAc,  dAy)  in  (30,  68}  which  may  not  be  physically  justifiable. 

5.3  Analysis  of  ESPRIT  under  Small  Array  Per¬ 
turbations 

Let  the  covariance  matrix  be  R  =  E{zz^},  where  z  =  and  let  its  eigen- 

y 

decomposition  be  R  =  E,A,E^  +  a'^EnE^.  The  ESPRIT  algorithm  solves  for  ® 
from  Er'®  =  Ey  and  estimates  the  DOA  angles  from  its  eigenvalues.  Here  Ex  and  Ey 
are  the  signal  subspace  eigenvector  matrices  for  the  X  and  Y  subarrays,  respectively, 
[50],  see  Chapter  4. 

Under  the  general  model  of  Eq.  (5.4),  the  application  of  the  ESPRIT  algorithm 
will  result  in  perturbation  of  the  eigenvalues  of  9  eis  the  next  theorem  shows. 

Theorem  1  Given  the  general  formulation  of  the  perturbed  array  signal  processing 
problem  of  Eq.  (5.4),  the  perturbed  k~th  eigenvalue  of^  can  be  written  as 

ik  =  4*  +  where  =  4*bfc(9Ay  -  9Ax)efc.  (5.11) 

Here,  =  exp  j— =  eJ(rA)'‘’  and  ejt  is  defined  as  the  n  x  1  unit  vector 
with  unity  value  at  the  k-th  element  position  and  zero  elsewhere. 

Proof  The  signal  subspace  eigenvector  matrix  E,  spans  the  perturbed  steering  ma- 


r  ^  gx  )  1 

trix  .  Thus,  there  exists  a  non-sineular  matrix  T  such  that: 

(FA  +  dAy)^ 


E,= 


(FA  +  ^A^) 
(FA  +  ^Aj)*^ 


(5.12) 


Solving  for  the  least  squares  solution  of  $  and  using  a  first  order  approximation  for 
the  pseudo-inverse  of  (FA  +  dAx),  [19],  see  also  Section  3.5.2,  yields 


9  =  T-'  {l  +  (rA)+(8A,  -  dA.)}  *T. 


(5.13) 


Since  the  eigenvalues  of  a  matrix  are  invariant  under  a  similarity  transformation,  the 
eigenvalues  of  4^  are  the  same  as  the  eigenvalues  of 


{l-h(FA)+(aA„-9A,)}^ 


(5.14) 


Using  first  order  perturbation  approximation,  [15]  and  see  also  Section  3.5.1,  for  the 
perturbed  fc-th  eigenvalue  yields 


6  =  6  +  6enrA)+(5A,  -  dAx)ek^ 


(5.15) 


The  proof  is  done. 


Having  found  the  perturbed  eigenvalues  of  4',  the  next  step  in  ESPRIT  is  to  solve 
for  the  estimated  DOA  angles.  The  next  lemma  gives  a  convenient  form  for  evaluating 
the  mean  squared  error  of  the  fc-th  estimated  DOA  angle. 


Lemma  1  Given  the  perturbed  eigenvalue  ik  =  4-  the  mean  squared  error  of 


the  corresponding  estimate  of  the  k-lh  DOA  is 


(5.16) 
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Proof:  Note  that  ~  and  that  ^k  =  ^k  +  ^^k  =|  It  1  where 

=  ^t  +  A^t  Using  first  order  approximation  and  sin{6k)  —  sin{dk)  +  A0fccos(0it) 
we  get 

A0fc|  (5.17) 

where  Ar*  denotes  the  perturbation  in  the  magnitude  of  the  estimated  eigenvalue  ^jt. 
Discarding  all  term  higher  order  terms  except  for  first  order,  we  get 

=  (5.18) 

Solving  for  |  A^t  yields  the  desired  result. 

Using  Lemma  1  and  Theorem  1  we  get  the  following  corollary. 

Corollary  1  The  mean-square~error  of  the  ESPRIT  doa  estimates  under  model  er¬ 
rors: 

^‘1  I’’  =  5  (2s)'  ■  (5-19> 

Let  Dt  =  Diag{Bk},  a  diagonal  matrix  formed  from  the  elements  o/a^  and  the  matrix 
terms  Mot;  Mit  be  specified  as  follows. 

Errors  in  Sensor  Gains  and  Phases 

Mot  =  Dt£{ vtvf  and  Mu-  =  DtE{ vtv[ }D[.  (5.20) 

where 

vt  =  [e^^MAdi  -  AQi)-fia,e^'^>(A^i  -  Aij^i).. 

•..e^-^'"(Ad;„  -  Aq,„)  +;a;;,e^*-(A<^;;,  -  A<i;„)]^.  (5.21) 

Errors  in  Sensor  Locations 


6  +  A^t  =  (l+Art)e-J^*‘"^» 


r2Kd 

j——cosek 


Mot  =  PDfcEfvtvf  }(rDt)"  and  Mu  =  rDtE{vtv[}(rDt)^. 


(5.22) 
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where 


Vfc  =  -Jy  -  ^H\}sinOk  +  (^i'l  -  )cos0it.. 

■■■{Afim  -  AfimjstnOk  +  {Ai>m  -  Ai/m}cosdk]^ .  (5.23) 


Errors  in  Sensor  Pair  Alignments 

Mofe  and  Mu  are  as  defined  tn  Eq.  (5.22)  where 

.  27tdcosdk  I 


Vfc  =  +j- 


(5.24) 


Proof:  Since  the  proofs  for  the  different  error  scenarios  follow  the  same  line  of 
argument  we  shall  present  only  the  proof  for  the  sensor  gain  and  phase  error  model. 
Note  however,  that 

dAr  =  AFxA  and  dAy  =  AP^A  (5.25) 

Using  Theorem  1,  the  perturbation  in  the  k-th  eigenvalue  is 

A(k  =  6bl.(Ar,  -  Ar,)Aefc 

=  6bl.(Ar,  -  Ar,)a*  (5.26) 

Note  however  that  defining  and  Djt  as  in  the  corollary  and  Eq.(5.21),  we  can 
re-write  the  above  as 

=  ^fcb5.DfcVfc  (5.27) 

Apply  Lemma  1  and  the  proof  is  done. 

Remark:  In  general,  the  matrix  terms  Mo*.-,  Mu-  are  not  diagonal.  However  in 

the  case  where  the  X  and  Y  subarrays  do  not  overlap  and  the  sensor  errors  (i.e., 
gain/phase,  location  or  alignment  errors)  have  zero  means  and  are  uncorrelated,  these 
terms  are  found  to  be  diagonal  and  more  explicit  expressions  for  Mofc.Mu  can  be 
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derived.  The  next  corollary  follows  from  Corollary  1  under  the  non-overlapping, 
uncorrelated  error  assumption. 


Corollary  2  Suppose  that  the  subairays  Ho  not  overlap  and  that 

1.  The  errors  in  sensor  gains  and  phases  are  :ero  mean,  uncorrelated  with  variance 
crl  and  aj,  respectively.  Then 

Moi  =  2Diag{eTl  +  •  •  • ,  (5-28) 

and 

Mu  =  2Diag{e^^*'  (a*  -  (tIq]),  •  •  • .  (5.29) 

2.  The  errors  in  sensor  locations  are  zero  mean,  uncorrelated  with  variance  and 
al  respectively.  Then 


(t)  rr" 

(5.30) 

*  (t)  crlcos'^9k)  TT^. 

(5.31) 

3.  The  angle  errors  in  sensor  pair  alignments  are  zero  mean,  uncorrelated  with  vari¬ 


ance  (t|.  Then 


I  2Trdcos0k  \ 

Moi  =  I - j-Z  1  alTT^ 


(5.32) 


(5.33) 


5.4  Discussion  of  Important  Special  Cases 


In  the  subsequent  discussion,  it  will  be  assumed  that  the  nominal  sensor  gains  and 
phases  are  unity  and  zero,  respectively.  The  first  part  of  the  discussion  deals  with 


61 


the  one  and  tow-source  special  cases  under  an  arbitrary  geometry  array  (the  nominal 
array  geometry  must  still  fulfill  ESPRIT  invariance  requirement  of  course).  The 
second  part  of  the  discussion  is  concerned  with  an  arbitrary  number  of  sources  under 
the  assumption  of  large  m  (number  of  sensors)  and  ULA  structure  for  the  array. 

5.4.1  Arbitrary  Array  Geometry  -  One  and  Two  Source 
Case 


For  non-overlapping  subarrays,  simple  closed  form  solutions  for  the  one  source  case 
and  approximate  solutions  for  the  two  source  case  can  be  found  for  the  errors  as  in 
Eqs.  (5.28)-(5.33)  and  Eq.  (5.19).  Note  that  for  the  one  source  ceise,  it  is  easily  shown 
that 

I  bi  1'*=  1  and  i?e{b[D?b,}  =  — .  (5.34) 

771  m 

Proof:  The  pseudoinverse  of  the  one-source  steering  matrix  aj  can  easily  be  shown 
to  be  equal  to  Thus  the  quantities  in  Eq.(5.34)  are  easily  found. 

For  the  two  source  case,  the  MSE  is  expressed  approximately  in  a  simple  form 
when  the  number  of  sensors  in  a  subarray  m  is  large  enough  such  that  |  (<C  m. 

I  bfc  |^~  —  and  /?e{bfD^.bfc}  ~  .  k  =  1,2.  (5.35) 

m  m 


proof:  The  evaluation  of  the  pseudoinverse  of  a  two-source  steering  matrix  A  is 
straightforward  albeit  tedious.  After  some  manipulation  it  can  be  shown  that 


aja; 

m 


a^ai 

m 


(5.36) 

(5.37) 

(5.38) 


Next  assuming  that  |  |<C  m  the  quantities  of  Eq.(5.35)  are  found. 
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We  now  turn  to  the  special  case  of  a  uniform  linear  array  (ULA)  and  where  the  total 
number  of  sensors  tuq  is  assumed  to  be  large.  There  are  two  cases  for  the  selection 
of  the  ESPRIT  subarrays  discussed  here.  The  subarrays  are  chosen  from  the  ULA  to 
be  either  non-overlapping  (i.e.  ,  m  =  mo/2,  the  so-called  interleaved  array.  [35])  or  to 
have  maximum  overlapping  (i.e.  ,  m  =  mo  —  1),  with  the  distance  spacing  between 
sensors  in  a  sensor  pair  kept  at  d. 
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Before  continuing  on  with  the  evaluation  of  the  various  MSE  expressions,  we  first 
recall  from  Appendix  G  of  [62]  that  for  a  ULA  with  m  sensors  that 


lim  -^A^A  =  I 

(5.42) 

m— *00  Yfl 

Thus  this  implies  that 

b^  =  e[.(A^A)~‘A^  ~  -^aj[ 

(5.43) 

m 

for  large  m.  Also  note  that 

II 

Q 

(5.44) 

Having  found  these  useful  quantities  the  examination  of  the  MSE  for  the  non- 
overlapping  and  maximally  overlapping  subarrays  begins. 

Non-Overlapping  Subarrays 

The  non-overlapping  case  involves  the  evaluation  of  the  expressions  in  Corollary 
2  and  the  use  of  Eqs.  (5.43),  (5.44). 

For  the  case  of  zero  mean,  uncorrelated  errors  in  sensor  gains  and  phases  it  can 
be  shown  that 


b^Mofcb;  ~  +  (yl)  and  /?e{b[Mu-bfc}  - - ^^-2  -2' 


Hence,  applying  Theorem  1,  we  find 

E{\Aek  n  ~  ~  (2;rdcos<? J  ^■= 

Proof:  The  Mofe,Mifc  terms  are  found  from  Corollary  2  and  yields 


(5.45) 


(5.46) 


Mok  =  2{al  +  <71)1  and  Mu-  =  2((t^  -  «t^)I  (5.47) 


Using  Eqs. (5.43)  and  (5.44)  yields  Eq.(5.45). 


Similarly,  for  the  other  two  error  scenarios  we  get 
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Sensor  Location  Errors 

For  the  case  of  zero  mean  uncorrelated  errors  in  sensor  locations 

£{]  P)  ^  +  (5.48) 

mo  vdcosPfc/ 


Sensor  Pair  Alignment  Errors 


For  the  case  of  zero  mean  uncorrelated  errors  in  sensor  pair  alignments 


E{\A6k 


(5.49) 


Subarrays  with  Maximum-Overlapping 

The  case  with  maximum  overlapping  of  subarrays  is  handled  in  a  similar  manner. 
Note  however  that  the  evaluation  of  the  matrix  terms  Mo*,Mife  is  different  from  the 
Non-overlapping  case. 

Sensor  gain  and  phase  errors 


For  zero  mean,  uncorrelated  errors  in  sensor  gains  and  phases,  it  can  be  shown, 
when  the  subarrays  are  selected  such  that  they  have  maximum  overlapping  (I.e.  , 
m  =  mo  —  1),  that 


b^Mofcb;;  ~  4-  (tI)  and  Re{hlMu^k}  ^  -^{<^1  -  <^5),  •  (5.50) 

TTlr 


mf, 


Hence, 


E{\  n  ~ 


ml  \2irdcos0kj 


al,  k  = 


(5.51) 


Proof:  Consider  the  evaluation  of  Mo*.  It  requires  the  evaluation  of  E{vfcVfc}  with 
V*  as  defined  in  Eqs.(5.20),(5.21).  With  maximally  overlapping  subarrays  we  get 


2  -1  O 
-1  2  -1 

-1 

O  -12 


{<^1  +  <^l) 


(5.52) 
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Also, 

■  2  -1  O' 

E{-Vky[}  =  ^  ~  ./  (5-53) 

.0  -1  2  . 

Thus,  from  Eqs.(5.20)  and  (5.43)  we  have 

‘2-1  O' 

biM„»b;  ==  ^(1  ■  ■  ■  11  (<t’ -  ^5)(1  ■  ■  11'  (5.54) 

I  O  -1  2, 

This  is  equal  to  the  expression  for  b|.Moitb][.  found  in  Eq.(5.50).  b^Mubit  is  deter¬ 
mined  in  the  same  manner.  The  proof  is  done. 

Remarks:  Thus,  when  the  subarrays  are  chosen  with  maximum  overlappping,  the 
MSE  of  ESPRIT  estimates  are  smaller  by  a  factor  of  2mo  than  when  the  subarrays 
are  non-overlapping.  This  is  intuitively  satisfying  ^ls  in  this  case,  the  subarrays  have 
a  larger  aperture  than  in  the  non-overlapping  subarrays  case.  Again,  the  MSE  are 
affected  only  by  the  phase  errors. 

The  expressions  for  the  other  error  scenario  is  analogously 

Se  isor  location  errors 

For  uncorrelated  errors  in  sensor  locations  and  with  maximum  overlapping  subar¬ 
rays 

£{|  p}  ^  (5-55) 

Comparison  with  MUSIC  under  sensor  gain  and  phase  errors 

The  MSE  expressions  for  MUSIC  under  random  model  errors  can  be  found  in 
[67].  As  before,  by  arguments  similar  to  those  found  in  [62],  for  the  case  of  zero 
mean,  uncorrelated  errors  in  the  sensor  gains  and  phases,  the  MUSIC  MSE  is 
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Remarks:  It  can  be  seen  that  the  MUSIC  MSE  is  affected  by  both  errors  in  gains 
and  phaises.  In  contrast  ESPRIT  is  affected  only  by  phase  errors.  Furthermore,  from 
Eqs.  (5.46),  (5.51)  and  (5.56)  it  is  clear  that  the  MUSIC  MSE  is  smaller  when  the 
number  of  sensors  mo  is  large.  The  fact  that  the  MSE  of  MUSIC  is  less  than  that  for 
ESPRIT  for  ljurge  mo  is  also  consistent  with  the  results  in  [63].  Similar  observations 
for  ESPRIT  and  MUSIC  can  be  made  for  the  other  model  error  scenarios. 

5.5  Extension  of  Results  to  DOA-dependent  Sen¬ 
sor  Gain  and  Phase  Errors 

In  the  previous  discussions,  it  has  been  assumed  that  the  sensor  gain  and  phases  errors 
are  independent  of  the  DOA  angles  of  the  source  signals.  Surprisingly  all  the  previous 
results  concerning  the  MSE  errors  for  this  error  model  still  holds  when  this  model 
is  generalized  to  the  case  where  the  errors  are  dependent  on  the  DOA  angles  (i.e., 
the  errors  are  different  for  different  angles).  Specifically,  consider  the  case  where  the 
errors  in  gains  and  phases  are  different  for  different  DOAs.  Then  for  the  perturbed 
model  in  the  general  model  of  Eq.  (5.4),  we  write 

aA,  =:[Ar,ai  l-.-Arnan]  (5.57) 

and 

aA,  =  [Af,a,  i-.-Arnan]  (5.58) 

where  the  terms  AF/t,  AF^  are  the  first  order  perturbation  terms  for  the  k-th  DOA 
angle  of  the  X  and  Y  sub-arrays  respectively. 

Thus  Theorem  1  can  still  be  applied  using  these  generalisations  for  the  sensor  gain 
and  phase  errors  thereby  giving  the  k-th  perturbed  eigenvalue  from  Theorem  1  to  be 


Aa  =  6bUArt- Ar.)a. 


(5.59) 
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From  this  equation  we  can  see  that  it  is  in  the  same  form  as  the  DOA  independent 
error  model  and  hence  all  the  previous  discussions  regarding  the  sensor  gain  and  phase 
error  model  would  still  hold. 

5.6  Simulations 

Two  uncorreiated  sources  are  given  at  DOA  angles  —20®  and  —10°  both  with  unit 
power.  The  X  and  Y  subarrays  are  taken  from  a  uniform  linear  array  with  distance 
spacing  between  the  sensors  in  sensor-pairs  taken  to  be  d  =  The  nominal  sensor 
gains  and  phaises  are  unity  and  zero  respectively.  Zero  mean,  uncorrelated  errors  in  the 
sensor  gains  and  phases  are  introduced  with  standard  deviations  of  10~^  and  4  degrees, 
respectively.  The  total  number  of  sensors  mo  is  varied  from  6  to  20  and  a  hundred 
trial  runs  are  used  to  average  the  MSE  of  the  TLS-ESPRIT  (see  [50])  DOA  estimate 
for  both  the  non-overlapping  and  maximum  overlapping  case.  Their  respective  Root- 
MSE  are  plotted  in  Figs.5.1,  5.2  where  the  solid  line  depicts  the  approximate  Root- 
MSE  for  the  non-overlapping  case  as  computed  from  Eq.  (5.46),  the  dashed  line  gives 
the  approximate  Root-MSE  for  the  maximum-overlapping  case  as  computed  from  Eq. 
(5.51)  while  the  symbols  and  ’-f  ’  denote  the  simulation  Root-MSE  of  ESPRIT  for 
the  non-overlapping  and  maximum-overlapping  cases,  respectively.  The  approximate 
Root-MSE  for  MUSIC  as  computed  from  Eq.  (5.56)  shown  by  a  dotted-line  and  the 
simulation  Root-MSE  for  MUSIC  denoted  by  the  symbol  'o’  are  plotted  in  these 
figures. 


5.7  Conclusions 

The  MSE  expressions  of  the  ESPRIT  DOA  estimates  under  three  random  model  error 
scenarios  are  derived  assuming  that  the  covariance  matrix  is  known.  Solutions  for  the 
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Rooc-MSE  of  -20  deg.  OOA  estimate  vs.  No.  sensors 


Figure  5.1;  Estimate  of  DO  A  at  -20  deg, 

ESPRIT  MSE  with  an  arbitrary  ESPRIT  array  geometry  and  one  and  two  sources 
are  given.  Approximate  solutions  for  a  uniform  linear  array  (ULA)  and  arbitrary 
number  of  sources  with  a  large  number  of  sensors  are  also  given.  Solutions  for  both 
cases  suggest  that  ESPRIT  is  affected  only  by  phase  errors.  Furthermore,  it  is  shown 
that  ESPRIT  with  maximum-overlapping  subarrays  will  exhibit  lower  MSE  (which 
decreases  at  the  rate  ■^)  than  ESPRIT  with  non-overlapping  subarrays  (which  de¬ 
creases  at  the  rate  ^).  When  compared  to  MUSIC  estimates  (whose  MSE  decreases 
at  the  rate  ^),  it  is  found  that  MUSIC  estimates  will  generally  give  a  lower  MSE 
than  ESPRIT.  The  simulation  MSE  are  well  in  accordance  with  the  analytical  results. 
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A  SIGNAL  SUB  SPACE 
APPROACH  TO  MODEL 
ERRORS 


6.1  Introduction 

E/igenstructure  methods  such  as  MUSIC  [53]  and  ESPRIT  [50]  have  gained  consid¬ 
erable  interest  in  recent  years  due  to  their  superresolution  capabilities  in  resolving 
closely  spaced  sources.  These  methods  are  hampered  however  by  the  requirement  of 
a  known  array  manifold  (i.e.  ,  the  collection  of  array  signal  vector  for  all  possible 
direction  of  arrival  (DOA)  angles)  as  in  MUSIC  or  that  the  sensor  subarrays  as  in 
ESPRIT  be  matched. 

It  has  been  observed  that  in  practice,  even  after  hardware  calibration,  sensor 
arrays  still  suffer  from  gain  and  phase  errors,  location  errors,  etc.  ,  [40].  The  effects 
of  random  array  model  errors  on  the  DOA  estimates,  in  terms  of  their  mean-square- 
error  (MSE)  and  sensitivity,  of  algorithms  such  as  ESPRIT  and  MUSIC  have  been  a 
subject  of  many  recent  studies  (see  e.g.  [13]- [60]).  These  studies  have  shown  that  the 
mean-square-error  (MSE)  of  their  DOA  estimates  will  increase  due  to  these  errors. 
As  a  result,  several  methods  have  been  proposed  for  DOA  estimation  under  various 
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array  model  errors.  For  the  case  of  DOA  estimation  with  unknown  sensor  gains  and 
phases,  a  method  applicable  only  to  uniform  linear  arrays  is  proposed  in  [38].  In 
[79,  39]  methods  based  on  the  orthogonality  of  the  noise  subspace  to  the  steering 
vector  are  proposed.  This  algorithm  involves  the  computation  of  the  inverse  of  a 
matrix  that  is  close  to  being  singular  when  the  SNR  (Signal-to-Noise  Ratio)  is  high 
and  when  the  estimates  of  the  DOA  angles  are  close  to  the  true  ones,  as  will  be 
shown  in  Section  III  of  this  paper.  Calibration  methods  using  known  DOA  angles 
and  experimental  results  on  an  ultrasonic  sensor  array  testbed  is  presented  in  [40].  An 
approach  requiring  uncorrelated  sources  is  found  in  [6]  while  a  maximum-likelihood 
based  algorithm  which  requires  known  field  covariance  (i.e.  .  the  covariance  matrix 
with  no  sensor  gain  and  phase  errors)  is  proposed  in  [14]. 

Studies  of  the  Cramer- Rao  lower  bounds  of  the  DOA  estimates  in  the  presence 
of  random  sensor  locations  errors  and  an  algorithm  for  DOA  estimation  under  these 
errors  can  be  found  in  [47],  [48].  The  algorithm  in  [48]  requires  that  sources  do  not 
overlap  either  in  time  or  in  frequency.  An  algorithm  for  estimating  unknown  sensor 
locations  is  proposed  in  [31]  but  requires  that  the  DOA  angles  be  known.  An  iterative 
maximum- likelihood  algorithm  was  proposed  in  [SO]. 

In  this  paper,  a  constraint  on  the  set  of  sensor  gains  and  phases  and  DOA  angles 
consistent  with  the  signal  subspace  is  studied.  This  constraint  is  shown  to  yield 
a  simple  way  of  computing  the  sensor  gains  and  phases  for  a  given  set  of  DOA 
etngles.  The  issue  of  uniqueness  of  the  set  of  sensor  gains  and  phases  and  DOA 
angles  is  examined.  It  is  shown  that  when  the  DOA  angles  are  known,  the  estimates 
of  the  sensor  gains  and  phases  computed  using  the  constraint  are  asymptotically 
unbiased.  When  the  DOA  angles  are  unknown,  an  iterative  procedure  incorporating 
the  constraint  is  proposed.  The  iterative  procedure  is  shown  to  be  able  to  handle 
other  array  error  models,  notably,  errors  in  sensor  locations  and  errors  in  sensor-pair 
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alignments  of  ESPRIT-type  subarrays.  Simulation  results  are  presented  to  show  that 
the  iterative  procedure  outperforms  conventional  ESPRIT  under  array  model  errors. 

6.2  Definitions  and  Model  Formulation 

The  following  notational  convention  is  used  in  this  paper.  .All  matrix  quantities  will 
be  denoted  by  bold  faced  capital  letters  while  vector  quantities  will  be  denoted  by 
bold  faced  small  letters.  (.)^,  (.)^,  (.)“*  denote  matrix  transpose,  matrix  conjugate 
transpose  and  matrix  inverse  respectively.  The  determinant  of  a  matrix  is  denoted  as 
Det{.). 

The  vector  norm  used  here  is  the  Euclidean  vector  norm  denoted  as  j]  .  H  and  the 
matrix  2-norm  will  be  denoted  as  ||  .  ||2.  The  2-norm  of  a  matrix  M  is  defined  as  that 
induced  by  the  Euclidean  vector  norm,  [15j: 

11  M  Il2=  ^maXo  II  Mu  II  (6.1) 

The  notation  Diag{v}  is  taken  to  denote  a  diagonal  matrix  formed  from  the 
elements  in  the  column  vector  v. 

VVe  now  define  the  nominal  model  of  the  DO  A  estimation  problem. 

Let  m  be  the  number  of  sensors  in  the  array,  n  be  the  number  of  narrowband  far- 
field  sources  and  (x,-,  yi]  be  the  location  of  the  i  —  tu  sensor  in  the  array.  The  nominal 
sensor  gains  and  phases  are  assumed  to  be  equal  to  unity  and  zero,  respectively.  Then 
the  nominal  model  for  the  nar.  3wband  array  signal  processing  problem  can  be  stated 
as  follows; 

z—Xs  +  y  (6.2) 

where  z  denotes  the  observation  vector  from  the  sensor  array,  A  is  the  steering  matrix 
defined  by  the  set  of  distinct  DOA  angles  {<?*}  (the  DOA  angles  are  measured  with 
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reference  to  the  y-axis  of  the  sensor  array). 

A  =  I  a(t^,)  •  •  ■  I  .  (6.3) 

Here,  the  nominal  steering  vector  at  DOA  angle  0  as 

27r 

a(6?)  =  [I  exp{—j -^{xiStnd  +  yycosd)} ... 

A 

'2‘!r  x 

4- i/m-icos5)}]  (6.4) 

where  A  is  the  wavelength  of  the  narrowband  signals,  s  is  the  source  vector  and  x] 
denotes  the  additive  noise  at  the  sensor  array. 

The  following  assumptions  are  made  throughout  this  paper. 

(A.l)  The  array  of  m  sensors  is  unambiguous.  [34],  i.e,  ,  it  satisfies  the  prop¬ 
erty  that  for  any  collection  oi  k  <  m  distinct  DOA  angles,  the  matrix 

[a{^i),  •  •  • ,  a(dfc)]  has  full  column  rank,  where  m  >  n,  and  n  is  known. 

(A.2)  The  source  covariance  matrix  Q  =  £^{ss^}  is  positive  definite. 

(A. 3)  s  and  p  are  mutually  uncorrelated  and  are  stationary  processes. 

(A.4)  T]  is  additive  white  Gaussian  with  covariance  £’{7777^}  =  ct^I. 

From  the  observation  vector  z  the  nominal  covariance  matrix  can  be  expressed  as 

Ro  =  E{zz^]  =  AQA^  +  aH.  (6.5) 

6.3  Array  Processing  Under  Unknown  Sensor  Gains 
and  Phases 

III.l  Coir.putation  of  Sensor  Gains  and  Phases 

For  unknown  sensor  gains  and  phases  the  array  model  would  be  written  as; 


z  =  FAs  -f  7 


(6.6) 


with  r  being  the  dUgonal  matrix  of  unknown  i  DOA  angle  independent)  sensor  gains 
and  phases  where 


r  =  Diag{  1,  aie-'"*’* .  •  -  • ,  } 

with  a,  and  <^,  denoting  the  t-th  sensor  gain  and  phase,  respectively.  We  assume  that 
all  sensor  gains  are  non-zero.  Without  loss  of  generality,  the  first  diagonal  element  in 
r  is  assumed  to  be  unity. 

Computing  the  covariance  matrix  of  z  yields 

R=  £;{2zf}  =  (rA)Q(rA)^  ^cr-I.  (6.7) 

Taking  the  eigenvector-eigenvalue  decomposition  of  R  gives 

R  =  E,A,El  -I-  a^E^El  (6-8) 

where  E,  denotes  the  eigenvector  matrix  of  R  associated  with  the  signal  subspace 
spanned  by  FA,  A,  is  the  eigenvalue  matrix  associated  with  it  and  En  denotes  the 
eigenvector  matrix  associated  with  the  noise  or  null  subspace  with  eigenvalue  cr^. 
Thus,  the  solution  set  of  sensor  gains  and  phases  and  DOA  angles  are  constrained 
such  that 

EsE^FA  =  FA,  (6.9) 

The  constraint  of  Eq.(6.9)  is  equivalent  to 

B,E*ra(«,)  =  ra(«,),i=  l,2,.-,n.  (6.10) 

By  defining  D,  =  Diag{a(^,)},  a  diagonal  matrix  composed  of  the  components  of  the 
steering  vector  a(5,),  and  v  as  the  m  x  1  column  vector  formed  from  the  diagonal 
components  of  F,  i.e.  ,  v  =  [1  •••  ]^,  we  can  rewrite  Eq.(6.10)  as 


E,EjDiV  =  D.v,  r  =  1, 2,  •  •  • ,  n. 


(6.11) 


Since  D,  is  a  diagonal  matrix  composed  eniireiy  of  the  steering  vector  at^;)  (see 
Eq.(6.4)),  it  is  also  unitary,  i.e.  D|D,  =  D,D|  =  I  and  therefore 

DiEsEp.v  =  v.z  =  (6.12) 

In  other  words,  from  Eq.(6.12),  the  vector  of  the  sensor  gains  and  phases  must  be  an 
eigenvector  with  unit  eigenvalue  of  the  matrices  W,  as  defined  by 

Wi  =  DfEsE^D,.  (6.13) 

We  note  that  W,  has  n  unity  eigenvalues  and  m  —  n  zero  eigenvalues.  Observe  that 
the  constraint  of  Eq.(6.9)  can  now  be  re-stated  as  the  set  of  n  eigenvector  constraints 
of  Eq.(6.12).  This  set  of  constraints  will  be  used  in  computing  the  sensor  gain  and 
phase  vector  v.  The  following  theorem  shows  how  these  n  eigenvector  constraints  can 
be  combined  into  a  single  eigenvector  constraint. 

Theorem  2  Let  D,-  and  W,  be  defined  by  Eqs.  (6.11)  and  (6.13),  respectively. 
Then  (V,  1)  is  an  eigen-pair  (eigenvector-eigenvalue  pair)  for  W,  =  DIEsEJD,, 
f  =  1,2,  ■  ■  ,n  if  and  only  if  (v,  n)  is  an  eigen-pair  for 

W  =  f;w.-  =  X^DjEsElD..  (6.14) 

1=1  «=i 

where  v  is  the  vector  defined  in  Eg.  (6.11).  Furthermore,  the  eigenvalue  n  is  the 
largest  eigenvalue  ofW. 

Proof:  It  is  e^y  to  see  that  if  v  is  an  eigenvector  with  unit  eigenvalue  for  the  matrices 
Wi  for  f  =  1,2,  •  •  • ,  n,  then  it  will  also  be  an  eigenvector  with  eigenvalue  equal  to  n 
for  W. 

So,  we  need  only  show  that  if  v  is  an  eigenvector  of  W  with  eigenvalue  n,  then  it 
is  also  an  eigenvector  of  W,  with  unit  eigenvalue. 
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Suppose  that  Wv  =  nv.  We  write  v  =:j  v  ii  v.  wnere  v  iias  unit  norm.  Then. 

•I 

11  Wv  il=  u  <  V  ii  W.v  (6.15) 

1=1 

Since  the  2-nortn  of  a  Hermitian.  semi- positive  definite  matrix  is  equal  to  its  largest 
eigenvalue,  we  get  ||  W,  1|2=  1.0.  From  Eq.fG.lo)  and  the  definition  of  the  matrix 
norm  as  in  Eq.(6.1),  we  conclude  that 

II  W.v|l=  1.0,2=  1, 2, (6.16) 
Eq.  (6.16)  is  sufficient  to  ensure  that  v  and  consequently  v  is  an  eigenvector  of  W, 
with  unit  eigenvalue. 

To  show  that  n  is  the  largest  eigenvalue  of  W,  we  need  to  show  that  |j  W  |j2=  n. 
Suppose  that  Wv  =  nv.  Then,  from  the  definition  of  the  2-norm  of  a  matrix  as  in 
Eq.(6.1),  we  get 

11  W  ||2>  n.  (6.17) 

However,  we  also  have 

11  W  |U<  x;  11  w,  Ib=  n.  (6.1S) 

1=1 

Therefore,  j|  W  ||2=  n. 

Remark:  Theorem  2  provides  a  way  to  compute  the  sensor  gains  and  phases  given 
the  DOA  angles,  utilizing  the  signal  subspace  constraint  of  Eq.  (6.9).  .As  such, 
calibration  of  the  array  for  sensor  gains  and  phases  can  be  done  for  known  DOA 
angles. 

The  following  corollary  is  an  immediate  consequence  of  the  previous  result. 

Corollary  3  The  eigenvector  v  found  in  Theorem  2  is  also  an  eigenvector  associated 
with  the  zero  eigenvalue  ofTS/l,  where 

M  =  f;DlE„Ej,D.. 

i=l 


(6.19) 


Remarks:  In  [79),  an  alternative  approach  was  taken  to  arrive  at  related  results. 
The  approach  of  [/9]  is  based  on  the  orthogonality  of  FA  to  the  noise  subspace  and 
a  solution  for  v  was  proposed  as  : 


M-‘eo 

ejM~‘eo‘ 


(6.20) 


where,  eo  =  [1  O---  0]^.  It  is  seen  from  the  previous  corollary  that,  at  the  true 
DOA  angles,  the  matrix  M  is  singular.  In  simulations,  it  has  been  observed  that 
the  estimate  of  this  matrix  tends  to  have  a  large  condition  number,  although  it  is 
usually  still  invertible  (due  to  noise,  errors  in  DOA  estimates,  etc.  ).  In  (39)  the 
orthogonality  of  FA  to  the  noise  subspace  is  also  exploited.  From  the  requirement 
that  E^ra(d<:)  =  0,k  =  1,2,  ..,n,  the  minimum  singular  vector  of  the  matrix  G  is 
proposed  as  the  the  gain/phase  vector  estimate,  where 


G  = 


E|,Di 


(6.21) 


EID„  J 

Note  that  the  evaluation  of  the  large  matrix  G  as  in  [39]  is  not  neccessary.  .Note  that 
the  right  singular  vectors  of  G  are  the  eigenvectors  of  G^G  and  therefore. 


G^G  =  Yi  D.^EnElD,  =  M  (6.22) 

i=l 

Thus  from  Corollary  3,  the  proposed  estimate  of  the  gain/phase  vector  of  [39]  and 
the  method  suggested  by  Theorem  1  are  equivalent. 

At  this  point  however,  it  is  not  immediately  apparent  that  the  computed  eigen¬ 
vector  V  is  unique.  This  issue  is  investigated  in  the  sequel. 

III.2  Uniqueness  Issues 


Given  the  signal  subspace  constraint  of  Eq.  iti.O).  there  are  three  questions  that 
should  be  addressed. 

First,  given  A,  is  the  corresponding  sensor  gain  and  phase  matrix  T  as  constrained 
by  Eq.  (6.9)  unique?  Second,  given  the  sensor  gain  and  phase  matrix  F.  is  the 
corresponding  steering  matrix  A  unique?  Third,  is  there  only  one  pair  of  sensor  gain 
and  phase  matrix  F  and  steering  matrix  A  which  are  consistent  with  the  constraint 
of  Eq.  (6.9)? 

We  examine  the  first  question  in  the  folowing  result  which  gives  a  sufficient  con¬ 
dition  for  the  affirmative. 

Theorem  3  If  m  ^  2n  —  I  and  if  all  size  n  subarrays  of  the  array  are  unambiguous, 
then  for  a  given  A  (as  defined  by  the  given  set  of  DO  A  angles),  the  set  of  senso'^  gains 
and  phases  as  constrained  by  Eq.(6.9)  is  unique. 

Proof  Appendix. 

Remarks 

1.  If  the  condition  in  Theorem  3  are  fulfilled,  the  matrix  W  of  Theorem  2  would 
yield  only  one  eigenvector  associated  with  the  eigenvalue  n. 

2.  The  condition  in  Theorem  3  that  all  size  n  sub-arrays  be  unambiguous  imposes 
certain  constraints  on  the  geometry  of  the  array.  It  is  also  equivalent  to  requiring 
that  all  n  X  n  sub-matrices  of  A  be  non-singular.  For  the  case  of  two  sources,  this 
condition  is  readily  met.  The  determinant  of  any  sub-matrix  of  A  can  then  be  written 
as: 

f  )  ~  -  a-jOfei),  i  7*^  (6-23) 

\  dki  ak2  ) 

The  quantity  in  Eq.(6.23)  will  be  zero  if  and  only  if:  a,iat2  =a,2ajn-  After  some  alge¬ 
braic  manipulation,  this  implies  that  0^  =  and  hence  would  violate  the  <issumption 
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that  A  has  full  column  rank. 

The  second  question  is:  given  the  sensor  gain  and  phase  matrix  F.  is  the  corre¬ 
sponding  set  of  DOA  angles  (i.e.  A)  as  constrained  by  the  signal  subspace  constraint 
unique?  The  answer  is  also  affirmative  with  the  condition  of  (A.l)  proving  to  be 
sufficient  in  this  case.  This  follows  directly  from  the  proof  for  Theorem  1  in  [76]. 

The  two  questions  addressed  so  far  only  guarantee  that  there  is  a  one-to-one 
correspondence  between  the  set  of  sensor  gains/ phases  and  the  set  of  DOA  angles  such 
that  the  constraint  is  fulfilled.  It  leaves  open  the  question  of  how  many  distinct  pairs 
of  sensor  gains/phases  and  DOA  angles  tliere  are.  This  leads  to  the  final  question,  i.e. 
,  is  the  set  of  sensor  gains/phases  and  DOA  angles  consistent  with  the  constraint  Eq. 
(6.9)  unique?  This  question  is  still  an  open  one.  We  can  however,  state  a  neccessary 
condition  which  must  be  fulfilled, 

Neccessary  condition:  For  uniqueness  of  the  pair  (F.A)  as  constrained  by  Eq.  (6.9), 
it  is  neccessary  that  for  all  diagonal  matrices  A  and  all  steering  matrices  Ai  ^  .A^, 
we  have  At  ^  AA^. 

Proof  Suppose  there  exists  a  diagonal  matrix  A  such  that  Ai  =  AA2.  Thus,  if  the 
pair  (F,  Ai)  satisfies  the  signal  subspace  constraint  Eq.  (6.9),  then  the  pair  (FA,  Aj) 
which  is  distinct  from  the  previous  pair  would  also  satisfy  the  constraint. 

Remarks: 

1.  It  can  be  easily  shown  that  a  linear  array  would  violate  the  above  neccessary 
condition.  In  (38)  and  in  [79]  it  has  been  shown  that  a  linear  array  would  not  yield 
unique  solutions  for  the  set  of  sensor  gains  and  phases  and  DOA  angles. 

III.3  Calibration  of  Sensor  Array 

Known  DOA  angles 
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If  the  DOA  angles  are  known  a-prton,  then  Theorem  1  can  be  applied  directly.  The¬ 
orem  2  provides  sufficient  conditions  for  uniqueness  of  the  estimate  of  the  sensor  gain 
and  phase  vector.  Therefore,  the  eigenvector  associated  with  the  largest  eigenvalue 
of  the  matrix  W  (as  defined  in  Theorem  1)  is  taken  to  be  the  estimated  sensor  gains 
and  phase  vector.  Suppjse  that  iV  snapshots  of  the  observation  vector  z  are  taken 
and  the  covariance  matrix  is  estimated  as: 


R=  ■^^z(*)z(i)^  (6.24) 

i=i 

The  estimate  of  the  sensor  gain  and  phase  vector  in  Theorem  1  can  be  shown  to 
be  asymptotically  unbiased.  Using  first  order  perturbation  analysis,  the  estimated 
signal  subspace  eigenvector  matrix  can  be  written  as  E,  =  E,  +  AE,,  where  AE,  is 
the  perturbation  in  the  signal  subspace  due  to  the  finite  number  of  snapshots  and 
noise  effects.  Thus,  we  get  the  estimated  matrix  W  =  W  +  AW,  where 

N 

AW  =  ^  Dt(AE,Ej  +  E,AEl  +  AE.AEDD,-  (6.25) 

isl 


Let  the  eigenvalues  of  W  be  denoted  by  where  g\  =  n  >  gj  Qm  and 

their  associated  unit  norm  eigenvectors  by  {v^.}.  Thus,  we  have  vi  =  where  v  = 
[1  ate-'*’  •  •  •  QTO-ie-'*”"'  ]^.  Then  the  estimate  of  the  sensor  gains  and  phases  vector 
is  written  as 


vi  +  Avi 

V  =  -= - 

e^(v,  +  Av,)eo 

where  v  =  [1  dje-'*’  •  •  ■  dm-ie-'*'"-’  p,  eo  =  [1  0  •  •  •  0]^  and 


(6.26) 


m 

Avt  =  f; 


fcs:2 


{vjAWv,} 

- p-Vfc 


(6.27) 


Approximation  to  first  order  yields 


V  =  v+  II  V  II  Av,-  II  V  (I  v(ejAv,eo) 


(6.28) 
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Thus,  see  [23]  for  evaluation  of  E{AE,AE|}  and  Ef  AE,}.  we  obtain 

E{v}  =  v  +  0(iV-')  (6.29) 

Hence,  when  the  DOA  angles  are  known,  the  estimates  of  the  sensor  gains  and  phaises 
as  computed  from  Theorem  i  are  asymptotically  unbiased. 

Unknown  DOA  angles 

In  the  case  where  both  the  set  of  sensor  gains  and  phases  and  the  DOA  angles  are  not 
known,  an  iterative  procedure  which  utilizes  the  constraint  Eq.  (6.9)  is  proposed.  The 
algorithm  starts  with  the  initial  DOA  estimates  found  by  assuming  nominal  sensor 
gains  and  phases  and  applying  the  eigenstructure  algorithm.  This  algorithm  can  be 
either  MUSIC  or  ESPRIT  or  any  other  eigenstructure-based  algorithm.  The  iterative 
procedure  is  summarized  in  the  following: 

Step  I,  Given  the  previous  DOA  angle  estimates,  compute  the  matrix  W  where  W  = 
DjEgElD,.  Here,  D,  =  Diag{a(0,)}.  Evaluate  the  eigenvector  v  associated  with 
its  largest  eigenvalue  and  use  it  as  the  updated  estimate  of  the  sensor  gains  and  phases. 

Step  II.  Apply  an  eigenstructure  method  to  estimate  the  DOA  angles,  using  the  up¬ 
dated  sensor  gains  and  phases  to  modify  the  method  appropriately  where  F  =Dia.g{v}. 

Step  III.  Check  for  convergence  of  the  estimates  of  the  sensor  gains  and  phases  and 
DOA  angles  according  to  some  convergence  criterion.  If  converged,  stop.  If  not,  go 
to  Step  I. 

In  the  above  procedure.  Step  II  specified  for  the  use  of  ESPRIT  can  be  stated  as 
follows.  The  matrix  FA  is  defined  as 


is  the  rotational  invariance  diagonal  matrix.  We  partition  the  matrix  E,  = 

Then  solve  for  the  matrix  ^  from  the  equation: 

E.v^'  =  r.vrp'EK.  (6.30) 

The  DOA  angle  estimates  are  then  computed  from  the  eigenvalues  of  ’S',  [50j.  In 
Section  V,  the  use  of  ESPRIT  in  the  iterative  procedure  will  be  called  M-ESPRIT 
(Modified-ESPRIT). 

6.4  Array  Processing  Under  Other  Error  Models 

Under  appropriate  assumptions,  other  model  error  scenarios  can  be  cast  in  the  same 
form  as  £q.  (6.6).  One  assumption  that  is  needed  is  that: 

(B.l)  The  magnitudes  of  the  DOA  angles  lie  within  a  “narrow”  enough  angular 
region  such  that  their  cosine  values  are  approximately  equal  to  each  other,  i.e.  , 
cos$k  ^  cosOi  =  const,  ior  k  ^  i. 

Errors  in  Sensor  Locations  (see  [60]):  Under  small  errors  in  sensor  locations, 
denoting  the  perturbed  location  of  the  ^'-th  sensor  as  (i^  +  +  Ayt),  we  write 

the  model  as: 

z  =  As  +  q  (6.31) 

where  A  =  [a(0i),  -  •  •  ,a(5n)]  'vith  a(5fc)  =  J'ka{0k)  and 


,  ,27rstndt„. 

I  -j - - - Dtag{Axu-  ■  ■ ,  Ai^} 

.  2vcosdk  •  f  ,  A  1 

-J - : - Diag{Ayi,  •  ■ ,  Ay„}. 


(6.32) 


If  the  sensor  location  errors  are  such  that  Axfc  <  Ayt,  i  e.  the  error  is  much  less 
along  the  x-axis  than  the  y-axis  (which  is  likely  to  be  the  situation  in  an  underwater 


towed  array,  [67])  and  using  assumption  (B.l),  we  can  approximate  the  matrix  by 
a  constant  matrix  F  where 

r  ~  I  -  j~Diag{Ay,,  •  •  ■ .  ^ym]xconst. 

Thus  the  model  can  be  written  as  in  Eq.  (6.6). 

Errors  in  Sensor-Pair  Alignments  (see  [60]):  In  ESPRIT,  two  subarrays,  denoted 
here  m  X-  and  Y-subarrays  are  configured  in  such  a  way  that  their  observation  vectors 
are  related  by  a  rotational  invariance  diagonal  matrix.  Each  sensor  in  the  X-subarray 
is  paired  with  a  sensor  in  the  Y-subarray  and  all  of  these  sensor  pairs  are  aligned 
such  that  their  normal  face  the  same  direction,  (see  Fig.  1).  A  source  of  error  would 
then  be  the  misalignments  of  these  sensor-pairs  such  that  they  are  rotated  about 

their  mid-points  by  small  angles  {^k}k~i . m*  Here  m  is  tno  number  of  sensors  in  a 

subarray.  We  write  the  model  as: 

z  =  As  +  p  (6.33) 


where 


B.v 


Where  $  is  the  rotational  invariance  matrix  defined  by 


$  =Diag{cxp(-j2^sm^i), ....  exp{-j^sin9n)} 


and 


(6.34) 


Here 


Bx  =  [rxib(0,),---,r.v„b(^„)l 
By  =  [Fnb(0, ),•••,  FynMM 


(6.35) 
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and 


fy»  =  I+j2;^Diag(^„ 

Here,  d  is  the  spacing  between  the  sensors  in  a  sensor-pair  and  h{6k)  is  the  nominal 
steering  vector  of  the  X-subarray  at  angle  Ok. 

Under  assumption  (B.l),  we  can  write  the  error  model  in  the  same  form  as  Eq. 
(6.6)  with 


Tx  ~  I  -  j3^Diag{/?,,  •  •  • .  x  const. 

Ty  ~  I-b  j^Diag{/3i,-  •  -,0,^}  X  const. 

Thus,  the  problem  of  estimating  DOA  angles  under  the  scenarios  of  errors  in 
sensor  locations  (where  errors  in  the  y-axis  dominates)  and  where  the  sensor-pairs 
in  ESPRIT  type  of  subarrays  are  mis-aligned  can  also  be  handled  by  the  proposed 
iterative  algorithm  also. 

Remark:  The  situation  where  more  than  one  type  of  errors  are  present  can  still  be 
handled  by  the  iterative  procedure  as  their  combined  error  model  would  still  be  of 
the  same  form  as  that  of  the  single  error  model. 

6.5  Simulation  Examples 

Example  1  (Errors  in  sensor  gains  and  phases) 

The  ESPRIT  subarrays  used  in  the  simulation  experiments  are  configured  as  in 
Fig.  4.3.  The  total  number  of  sensors  in  the  array  is  eight,  i.e.  m  =  8,  and  the  spacing 


between  sensor  pairs  is  taken  to  he  d  =  t-  Nominal  sensor  gams  and  phases  are  unity 
and  zero  respectively. 

There  are  three  uncorrelated  sources  at  DOA  angles  of  -30,  -f  o,  and  4-'20  degrees, 
all  at  the  same  SNR  level.  The  means  and  MSEs  are  averaged  over  two  hundred 
trials.  ESPRIT  and  M-ESPRIT  are  applied  with  fifty  iterations  of  M-ESPRIT  per 
trial. 

Fig.  2  shows  a  typical  plot  of  the  MSE  vs.  SNR  for  the  estimate  of  the  DOA  at  +20 
degrees.  The  solid  and  dashed-lines  denote  M-ESPRIT  and  ESPRIT,  respectively. 
Four  hundred  snapshots  are  taken  per  trial  and  errors  are  introduced  such  that  the 
gains  are  perturbed  by  up  to  8  percent  of  their  nominal  values  while  the  phases  are 
perturbed  by  up  to  9  degrees. 

Next,  the  SNR  is  fixed  at  20  dB  and  the  snapshot  size  (N)  is  varied.  The  sensor 
phases  are  perturbed  by  up  to  ±17  degrees  (with  zero  nominal  phase)  and  the  sensor 
gains  by  up  to  8  percent  of  the  nominal  gain  of  unity.  Fig.  3  shows  a  typical  plot  of 
the  MSE  for  the  estimate  of  the  DOA  angle  at  +20  degrees.  As  before,  the  solid  and 
dashed-lines  denote  M-ESPRIT  and  ESPRIT,  respectively. 

.M-ESPRIT  is  able  to  provide  better  DOA  angh  estimates  than  conventional  ES¬ 
PRIT  as  can  be  seen  from  Figs.  2  and  3.  Our  experiments  also  suggests  that  the 
DOA  angle  estimates  of  M-ESPRIT  show  smaller  bias  as  opposed  to  ESPRIT  esti¬ 
mates  which  demonstrate  significant  bias. 

Example  2  (Errors  in  sensor  gains,  phases  and  locations) 

There  are  three  uncorrelated  sources  at  DOA  angles  of  -20,  +5  and  +20  degrees  all 
at  SNR  of  30dB.  Five  hundred  snapshots  are  taken  per  trial  and  the  means  and  MSEs 
of  the  DOA  estimates  averaged  over  two  hundred  Monte  Carlo  trials  are  evaluated. 
The  sensor  locations  are  perturbed  in  the  i/-axis  by  up  to  O.lA  while  no  error  is 


induced  in  the  a:-axis.  and  the  sensor  gains  and  phases  are  perturbed  by  up  to  iO 
percent  and  ±  90  df^grees  respectively. 

The  M-ESPRIT  algorithm  (two  hundred  iterations  per  trial)  and  the  conventional 
ESPRIT  algorithm  are  applied  and  the  estimated  means  and  root-MSE  of  the  DOA 
estimates  are  tabulated  in  Table  1.  As  can  be  seen  from  Table  1,  the  iterative  al¬ 
gorithm  is  able  to  give  better  DOA  estimates  under  these  combined  errors  than  the 
conventional  ESPRIT  algorithm. 

6.6  Conclusions 

The  problem  of  array  signal  processing  under  unknown  sensor  gains  and  phases  is 
studied.  A  signal  subspace  constraint  is  used  to  yield  a  simple  way  of  computing  the 
sensor  gains  and  phases  from  a  given  set  of  DOA  angles.  The  question  of  uniqueness 
of  the  set  of  sensor  gains  and  phases  and  DOA  angles  is  discussed.  When  the  DOA 
angles  are  known  a-priori,  the  estimates  of  the  sensor  gains  and  phases  are  shown  to 
be  asymptotically  unbiased.  For  the  case  where  the  DOA  angles  are  also  unknown,  an 
iterative  procedure  is  proposed.  The  approach  in  the  paper  is  shown  to  be  applicable 
to  two  other  model  error  scenarios  under  appropriate  assumptions.  These  scenarios 
are:  unknown  sensor  locations  and  unknown  rotation  of  sensor  pairs  in  ESPPJT  type 
subarrays.  Simulation  experiments  show  that  the  iterative  algorithm  provides  better 
DOA  estimates  than  does  the  conventional  ESPRIT. 
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Proof  of  Theorem  3  We  argue  by  contradiction.  Suppose,  for  a  given  A.  that  there 
exists  two  sets  of  sensor  gains  and  phases  as  defined  by  the  matrices  Fi  and  r2  such 
that  Fi  r2  and  they  both  fulfill  Eq.(6.9). 

Then  EsEjFiA  =  FiA  and  E9E|r2A  =  F2A  This  implies  that  rank[A|AA]  = 
n,  where  A  =r7*r2  =  Diag{  1,  Ai,  •  •  •,  }.  Thus, 

ranA:[A  —  AA|AA]  =  n  (6.36) 

Writing  out  this  matrix  explicitly,  with  A  =  [0^*],  we  have 

1  •••  1 

Aiflu  •••  Aifli^n 

(6.37) 

Consider  the  (n  + 1)  x  (n  + 1)  sub-matrix  of  Eq.(6.37)  found  by  taking  the  block  from 
the  first  row  to  the  (n  -b  l)-th  row  and  the  first  column  to  the  (n  -b  l)-th  column.  We 
can  write  this  as: 


[A-AAIAA]  = 


0 

(1  —  Ai)aii 


0 


(1  —  Ai)ai. 


1  (1  -^m— 1  )^m— l.l  (1  "• 


submairix[A  —  AA|AA| 


0 

(1  —  Aj )aji 


0 

(1  Ai  )ni,n 
( 1  ~  A„)an,n 


1 

Aion 


(6.38) 


L  (1 

This  sub-matrix  is  singular  from  Eq.(6.36).  Taking  the  determinant  and  using  the 
fact  that  every  n  x  n  sub-matrix  of  A  is  non-singular,  we  conclude  that 


n(l-A.)  =  0  (6.39) 

«=i 

Hence,  from  Eq.(6.39),  for  at  least  one  t,  it  must  be  that  A,  =  l.O.  Using  the  same 
reasoning  repeatedly  as  in  eqs.(6.37)-(6.39),  we  can  show  that  at  least  (m  —  n)  of  the 
Aj  terms  are  equal  to  1.0. 


We  next  show  that  the  remaining  (n  —  i )  A,  terms  are  also  equal  to  i.O.  Without 
loss  of  generality  we  can  assume  that  the  lirst  (/n  —  n)  A,  terms  have  been  shown  to 
equal  1.0.  We  now  consider  the  matrix  {AjAAj  with  these  (m  —  n)  unit  values  of  A, 
plugged  into  the  matrix. 

Consider  the  (n  +  1)  x  (n  +  1)  sub-matrix  of  [AJAA]  found  by  taking  the  block 
from  the  first  to  the  (n  -f-  l)th  column  and  the  (m  —  2n  -f  l)</i  row  to  the  (m  —  n  -f  l)</i 
row.  Note  that  the  assumption  that  m  >  2n  —  1  ensures  that  there  are  enough  rows 
such  that  this  is  possible.  Then 

Um— 2n+l.I  '  '  U,n_2n+l,n 

su6matrta:[A|AA]  =  (6.40) 

®Tn— n,l  ■  ■  n.n 

.  n+1,1  n+l,n  Urn— n+1.1  . 

Subtracting  the  last  column  vector  from  the  first  column  vector  in  Eq.(6.40)  we  get 
the  matrix: 

0  <Jm-2n+l,2  "•  am-2n+l.n  0^-2n+l,l 

0  Utn— fi,2  ■  *  '  ^m—n.n  Un»-.n,l 

-^m— fi+l  )®Tn— n+l,l  *^Tn— n+1.2  n+l,n  '^m— n+1  Om— n+l.l 

Taking  the  determinant  of  Eq.(6.41)  which  is  equal  to  zero,  and  using  the  fact  that 
any  n  x  n  sub-matrix  of  A  has  to  be  non-singular,  we  get  the  result  that  either 
Om-n+i,!  =  0  or  =  1.0.  Since  all  the  elements  of  A  have  unit  magnitude  as 

can  be  seen  from  Eq.  (6.4)  implying  0,  we  conclude  that  =  1.0. 

Using  similar  reasoning,  we  can  show  that 

Aj  =  1.0,  i  =  m  —  n -f  1, ...,  m  —  1.  (6.42) 

Hence,  we  can  write  A  =  I  or  equivalently  Fi  =  r2.  This  contradicts  our  initial 
assumption  that  Fi  ^  r2.  O 


True  DOA 

Mean  (cleg.j 

-20 

-f20 

ESPRIT 

-14.279 

11.690 

43.372 

M-ESPRIT 

-18.949 

2.052 

21.290 

True  DOA 

Root-MSE  (deg.) 

-20 

+5 

-1-20 

ESPRIT 

5.721 

6.690 

mmm\ 

M-ESPRIT 

1.060 

2.976 

1.311 
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A  BLIND  IDENTIFICATION 
APPROACH  TO  SIGNAL 
SEPARATION 

7.1  Introduction 

In  recent  years,  there  has  been  interest  in  wideband  array  signal  processing  -  i.e., 
the  processing  of  wideband  signals  obtained  from  arrays  of  sensors.  This  is  of  interest 
in  such  diverse  areas  as  in  radar  and  sonar  processing,  where  parameters  of  inter* 
est  include  the  direction*of-arrival  (DOA)  of  multiple  signal  sources,  their  respective 
power,  signal  copy,  time  delay  estimation,  or  as  in  speech  processing  where  the  de¬ 
sired  speech  source  signals  are  extracted  from  corrupted  or  mixed  recordings  from 
microphone  arrays  (the  so-called  ’cocktail  party  problem’);  or  as  in  biomedical  signal 
processing  where  EMG  recordings  from  multiple  sensors  are  processed. 

In  practical  applications,  robustness  is  often  a  major  problem  for  many  modern 
array  signal  processing  techniques.  Most  algorithms  assume  certain  array  structure 
with  known  or  accurately  calibrated  characteristics  such  as  sensor  gains/phases  and 
locations.  In  reality,  errors  always  present  in  actual  physical  systems  and  as  a  result 
these  algorithms  may  not  perform  as  well  or  even  fail  completely  (60), [68,  69j. 
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Consider  the  teleconference  scenario  where  an  array  of  microphone  is  set  before 
an  audience.  The  speakers  may  be  distributed  near  or  far  away  from  the  array  and 
there  may  be  interference  sources  such  as  fan  hum.  music,  etc.  In  this  situation  it 
is  not  realistic  to  assume  that  the  speaker  locations  are  known  a-priori.  Note  also 
that  the  relative  powers  of  the  speakers  would  be  stronger  when  they  are  directly 
facing  the  array  and  not  as  strong  otherwise.  This  scenario  shows  that  the  source 
locations  and  sensor  gains  cannot  be  assumed  to  be  known.  Thus  there  is  a  need  for 
new  algorithms  or  approaches  that  is  robust  to  these  types  of  uncertziinties. 

7.1.1  A  New  Approach  and  Its  Features 

In  this  chapter,  a  new  method  based  on  the  so  called  blind  identification  is  proposed 
for  the  source  estimation  (signal  copy)  problem.  By  exploiting  signal  characteristics 
such  as  statistical  independencies,  the  proposed  approach  allows  maximum  uncertain¬ 
ties  associated  with  sensor  array  and  signal  propagation  environment.  The  key  feature 
of  the  proposed  algorithm  is  its  robustness  against  the  following  model  uncertainties; 

•  Unknown  Sensor  Gains:  Complete  calibration  of  sensors  may  be  costly.  Sensors 
may  have  different  gain  patterns  in  different  directions.  While  most  methods 
assume  known  sensor  gains,  the  proposed  method  is  not  only  insensitive  to  gain 
uncertainties,  it  also  provides  estimates  of  array  sensor  gains. 

•  Unknown  Combinations  of  Near-Field  and  Far-Field  Sources:  The  proposed 
method  can  deal  with  near-held  and  far-field  sources  simultaneously.  It  offers 
source  location  estimates  for  near-held  sources  and  direction  of  arrival  estimates 
for  the  far-held  sources.  Such  capability  is  potentially  important  in  applications 
such  as  robotic  vision. 
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•  Unknown  Combinations  of  Narrowband  and  Wideband  Sources:  The  sources 
may  be  a  combination  of  wideband  signals  such  as  voices  and  narrow-band 
signals  such  as  humming  noise  with  unknown  frequency  locations.  The  proposed 
approach  is  able  to  estimate  each  source  and  in  fact  identify  their  frequency- 
characteristics  {narrow-band  or  wideband). 

•  Unknown  Source  Spectral  Characteristics:  The  sources  can  be  overlapping  or 
non-overlapping  in  their  unknown  power  spectra.  The  proposed  approach  fully 
utilizes  the  entire  freqency  band  to  obtain  joint  estimates  of  the  signal  location 
and  channel  parameters. 

•  Unknown  number  of  signals:  The  proposed  method  does  not  assume  that  the 
number  of  signals  is  known.  A  new  clustering-based  technique  is  proposed  to 
determine  the  number  of  sources.  Such  a  method  is  crucial  when  the  number 
of  sources  can  not  be  determined  at  a  particular  frequency.  As  an  example, 
consider  the  case  where  there  are  four  sources,  two  of  which  have  overlapping 
spectra  disjoint  from  the  other  two.  Then  any  estimation  of  the  number  of 
sources  at  any  particular  frequency  would  yield  less  than  four  sources. 

The  proposed  method  will  be  demonstrated  both  for  wideband  scenarios  such  as 
separate  multiple  voice  sources  and  for  narrowband  plus  wideband  scenarios,  involving 
near-field  and  far-held  sources  and  different  array  uncertainties. 

7.1.2  Related  Works 

In  radar  and  sonar  applications,  the  directions-of-arrivals  of  multiple  sources  are  of 
much  interest.  A  variety  of  eigenstructure-based  methods  have  been  developed  for 
estimating  DOAs  of  these  sources  and  when  the  model  assumptions  of  these  methods 
are  satisfied,  they  are  capable  of  giving  good  estimation  performance,  see  eg.,  [53,  50]. 
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However,  in  the  presence  of  sensor  errors  or  model  uncertainties,  these  methods  would 
deteriorate  or  fail  altogether. 

Sensor  errors  in  the  actual  physical  array  system  can  arise  from  changes  in  the  envi¬ 
ronment  in  which  the  sensor  array  is  placed,  the  degradation  of  electronic  components, 
calibration  errors  which  result  in  uncertainties  in  the  sensor  gains  and  phases.  .As  a  re¬ 
sult  of  this,  several  eigenstructure- based  approaches  have  been  proposed  for  the  DOA 
estimation  problem  under  sensor  errors  for  narrowband  sources,  see  (81,  SO,  6].  These 
methods  are  usually  iterative  in  nature  and  involve  some  sort  of  multi-dimensional 
search  over  the  space  of  the  unknown  sensor  parameters  and  DOA  angles.  Conse¬ 
quently,  these  methods  are  computationally  intensive  and  may  suffer  from  convergence 
to  local  minimum  solutions. 

Microphone  array  processing  have  been  proposed  for  speech  de-reverberation,  en¬ 
hancement  and  noise  reduction,  co-channel  interference  rejection  and  extraction  of 
desired  speech  sources  from  noise  or  jammer  contaminated  observations,  see  (25,  36, 
32,  12,  16,  57,  22].  The  use  of  microphone  arrays  has  a  host  of  intriguing  and  excit¬ 
ing  applications  such  as  a  front  end  for  automatic  speech  recognition,  teleconferences 
and  hands-free  cellular  telephony  in  cars.  Most  of  these  approaches  are  based  on 
the  beamformer  concept  in  array  signal  processing,  where  the  individual  outputs  of 
the  array  are  delayed  and  summed  in  such  a  way  as  to  produce  a  narrow  beam  of 
sensitivity  towards  the  direction  of  the  desired  signal.  These  approaches  however  re¬ 
quire  model  assumptions  that  are  not  easily  satisfied  such  as  identical  microphones, 
ideal  point  sources  for  the  source  signals,  calibration  of  the  microphone  array,  known 
source  locations  or  in  the  far  field,  etc. 

This  chapter  is  organized  as  follows.  After  piesenting  the  problem  models  and 
assumptions  in  Section  2,  we  summarize  certain  results  in  blind  identification  and 
estimation  of  independent  sources  in  Section  3.  In  Section  4.  we  outline  the  proposed 
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approach  and  give  detailed  discussions.  We  present  simulations  separating  multiple 
voice  sources  to  demonstrate  our  approach.  We  conclude  with  some  comments  on 
some  variations  and  future  work. 

Throughout  this  chapter,  we  assume  the  following  notational  convention;  All 
matrix  quantities  will  be  denoted  by  upper  case  bold  faced  symbols  while  vector 
quantites  will  be  denoted  by  lower  case  bold  faced  symbols.  Also,  the  symbols 
(.)^,  (.)~^  and  {.)*  will  denote  the  transpose  of  (.),  the  complex  conjugate  transpose 
of  (.),  the  inverse  of  (.)  and  the  pseudo-inverse  of  (.)  respectively. 

7.2  Models  and  Assumptions 

The  problem  of  estimating  multiple  sources  is  illustrated  in  Figure  7.1.  From  the 


''th  Sensor 

^”*®ngure  7.1:  Estimation  of  Multiple  Sources 
signals  collected  by  m  sensors  our  objective  is  to  estimate  the  n  sources 


96 


{s,(<)}.  A  wideband  model  is  given  in  time  domain  by 

n 

-  r,fc)  +  i  =  1,2,  •  •  •  ,m.  (7.1) 

fc=i 

where 

x,(f):  received  signal  at  the  ith  sensor. 

m:  the  number  of  sensors  (known). 

a,*:  (unknown)  channel  coefficient  that  characterizes  attenuation  of  the  kih  source  at 
the  ith  sensor  as  the  results  the  sensors  and  sources  directionality,  sensitivity, 
gains,  geometrical  spreading,  etc.  Without  loss  of  generality,  the  first  sensor  is 
taken  as  a  reference  and  au  =  1. 

Sk'-  the  A:-th  (unknown)  source  signal. 

r,fc;  (unknown)  time  delay  at  the  i-th  sensor  of  the  k-th  source.  Again,  the  delays  are 
taken  as  the  relative  delays  with  respective  to  the  reference  sensor  and  Txk  =  0. 
The  delays  are  related  to  the  path  difference  given  by  nk  =  where  d{k  is  the 
path  difference  of  the  k-th  source  to  the  i-th  sensor  (relative  to  the  reference 
sensor)  and  c  is  the  speed  of  propagation.  See  Figure  (7.1).  The  time  delays 
are  due  to  the  spatial  distribution  of  the  sensors  and  source  locations  (i.e.,  the 
DO  A  angles  6  k  and  ranges  Rk).  Propagation  of  the  signals  through  the  channel 
can  be  far-field  or  near  field  or  any  combination  thereof. 

n:  the  number  of  sources  (unknown). 

The  corresponding  frequency  domain  description  of  the  above  model  is  obtained 

by  taking  Fourier  transform  of  the  time  domain  model 


x(tu)  =  A(u;)s(u;)  +  tj{u}) 


(7.2) 
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where,  x{w),  s(u;),  and  t/(u;)  are  the  Fourier  transform  of  x{^),s{f)  and  T}{t),  respec¬ 
tively,  and 

1.0  •••  1.0 

dm  i 

The  problem  then  is  to  estimate  the  gains  a,k  (the  channel  coefficients),  delays  Tik 
(or  DOA  angles  Ok  and  ranges  Rk)  and  the  individual  signals,  Sk{t)- 

Remarks:  Considerable  simplification  of  the  problem  can  be  achieved  if  the 
sources  are  all  narrowband.  However,  most  signals  of  interest  (such  as  acoustic  ones) 
are  not  narrowband  signals.  Consider  the  following  case  for  a  speech  signal.  For  a 
speech  signal  with  a  bandwidth  of  4kHz,  the  time  delays  would  have  to  be  much  less 
than  .25  milliseconds  to  meet  the  narrowband  assumption.  At  a  sound  velocity  (in 
air)  of  340  m/s,  this  implies  that  the  path  lenghts  of  the  signal  to  any  two  sensors 
cannot  differ  by  more  than  .085m.  Thus,  in  order  for  the  narrowband  assumption  to 
be  valid,  all  the  sensors  within  the  array  must  lie  within  a  distance  of  8.5cm  or  less 
from  each  other!  It  is  therefore  essential  to  consider  the  wideband  case. 

The  key  assumption  used  in  the  proposed  approach  is  signal  (statistical)  inde¬ 
pendency.  Although  such  an  assumption  appears  rather  restrictive,  it  holds  approxi¬ 
mately  among  signals  from  different  sources.  It  is  also  important  to  note  that,  from 
our  simulations  of  using  different  speech  and  music  signals,  our  algorithm  performs 
well  even  if  there  is  no  proof  that  the  signals  involved  are  statistically  independent. 
Specifically,  the  assumptions  adopted  in  this  chapter  is  given  as  follows. 

Basic  Model  Assumption 

Al:  5,(t)’s  are  zero-mean  and  mutually  independent. 

A2:  »7,(t)’s  are  zero-mean  Gaussian  noise  independent  of  s,(i)’s. 
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7.3  Proposed  Approach 

The  central  idea  of  the  proposed  approach  stems  from  the  idea  of  applying  blind 
identification  techniques  to  the  frequency  domain  model 

x{u)  =  A{u})siuj)  +  T]{u>)  (7.4) 

where  matrix  A(w)  contains  the  information  of  source  location  and  channel  charac¬ 
teristics  as  specified  by  (7.3).  The  source  can  then  be  extracted  from  the  observation 
by  applying  least-square  methods  to  (7.4).  .A.  schematic  diagram  of  the  proposed 
approach  is  shown  in  Figure  (7.2). 


x(t) 


Signal  and  Channel 

Signal  and  Channel 

FFT 

X(a») 

Parameter  Estimatioi 

Parameter  Estimation 

Based  on  a  Single 

Based  on  All 
Frequencies 

Frequency 

(Oustering  Algorilhn ) 

Figure  7.2:  Schematic 


Diagram  of  Proposed  Approach 


We  must,  however,  pay  special  attention  to  certain  subtleties  due  to  the  fact  that 
the  spectra  of  the  source  signals  are  unknown.  Because  each  source  may  or  may 
not  have  energy  at  a  particular  frequency,  not  all  the  source  locations  and  channel 


parameters  can  be  estimated.  1  he  challenge  is  to  obtain  estimates  at  each  frequency 
and  combine  them  to  form  joint  estimates.  We  shall  next  discuss  the  functions  of 
blocks  in  the  schematic  diagram. 

7.3.1  Blind  Identification  and  Estimation  of  Independent 
Sources 

The  so-called  blind  identification  and  estimation  deals  with  the  problem  of  estimat¬ 
ing  source  signals  when  the  transmission  channel  is  unknown  [10.  70].  This  seemingly 
ill-posed  problem  is  approached  by  exploiting  the  statistical  properties  of  the  source 
signals.  In  the  context  of  signal  estimation  problem  considered  in  this  chapter,  the 
statistical  independency  of  different  sources  enables  us  to  allow  various  model  uncer¬ 
tainties  as  outlined  in  the  Introduction. 

The  basic  identification  and  estimation  (blind)  problem  considered  in  this  chapter 
involves  a  linear  model 

x(t)  =  As(0  +  n(f)  (7.5) 

where  x  is  the  observation  random  vector,  matrix  A  corresponds  to  the  channel 
characteristics,  s{t)  is  the  source  random  vector,  and  n(f)  is  the  measurement  noise. 
The  objective  of  blind  identification  and  estimation  is  to  identify  and  estimate  both 
A  and  s{t)  from  x(f).  In  particulai,  we  deal  with  the  folowing  questions: 

Given  x{t)  generated  from  independent  source  s{t),  nonsingular  channel  parameter 
matrix  A,  and  additive  Gaussian  noise  n(<)i 
Ql:  to  what  extend  A  and  s{t]  can  be  determined? 

Q2:  how  to  estimate  A  and  s(f). 

These  questions  have  been  resolved  recently  [73],  [74],  and  their  answers  are  sum¬ 
marized  in  the  following. 
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Theorem  4  The  channel  parameter  matrix  A  uniquely  determined  by  the  obser¬ 
vation  statistics  up  to  post  multiplications  of  a  diagonal  matrix  and  a  permutation 
matrix  if  and  only  if  all  except  perhaps  one  source  are  non-Gaussian. 

An  important  implication  of  the  above  theorem  is  that  the  waveforms  of  the  inde¬ 
pendent  sources  can  be  recovered  from  the  observation  without  knowing  the  channel 
as  long  as  they  are  non-Gaussian.  Although  scales  and  orders  of  these  signals  can  not 
be  determined  without  further  information,  it  is  hardly  a  problem  in  practice  since 
the  waveforms  of  these  signals  are  obtained.  Such  kind  of  identification  is  referred  to 
as  waveform-preserving  blind  identification. 

There  are  different  methods  of  determining  the  channel  parameter  matrix  from 
the  statistics  of  the  observation.  Using  the  fourth-order  moments,  Cardoso  proposed 
the  fourth-order  blind  identification  algorithm  (FOBI)  [10).  The  FOBI  algorithm  has 
the  restriction  that  the  source  signals  must  have  distinct  kurtoses  and  be  noise-free. 
A  slight  extension  can  be  obtained  when  the  noise  is  spatially  white  or  the  noise 
statistics  are  known  [70],  [61].  Algorithms  using  cumulants  of  arbitrary-orders  was 
presented  in  [74],  These  algorithms  eliminates  the  restrictions  on  the  source  and 
applies  to  arbitrary  unknown  Gaussian  source. 

7.3.2  Source  Location  Estimation  using  a  Single  Frequency 

At  a  fixed  frequency  wq,  we  can  estimate,  by  a  applying  blind  identification  method, 
K  columns  of  matrix  A(u;)  that  correspond  the  K  signals  having  energy  at  wo,.  Con¬ 
sequently,  we  obtain  the  estimates  of  the  gain  factors  {a,^}  and  the  phases 
corresponding  to  these  K  signals.  The  key  question  addressed  here  is  how  to  infer 
the  locations  of  the  K  sources  from  the  estimates  }. 

There  are  two  ways  to  go  about  estimating  the  source  locations  from  the  estimates 


of  },  found  by  using  blind  identificaiion.  The  first  approach  applies  when  a 

certain  uniqueness  condition  imposed  on  the  sensor  and  source  holds.  In  [59j,  several 
closed  form  solutions  for  the  source  locations  given  the  relative  path  differences  is 
proposed  that  may  be  readily  applied  to  the  present  problem.  Note  however,  that  the 
relative  path  differences  are  determined  from  the  estimates  of  {e"-'"  )  which  limit 

the  path  difference  values  (to  avoid  phase  wrapping)  such  that 

I  -  1  max{|  d.-fc  1}  <  TT.  (7.6) 

c 

Unfortunately,  it  is  impossible  to  verify  the  above  condition  since  d.k’s  are  unknown, 
more  useful  (sufficient)  condition  can  be  obtained  easily  as 

1  -  I  r  <  T,  (7.7) 

c 

where  r  is  the  maximum  radial  distance  of  the  sensor  array.  The  above  constraint 
can  of  course  always  be  satisfied  if  only  the  low  frequencies  are  used.  The  advantage 
of  this  method  is  its  simplicity  while  the  disadvantage  is  that  the  entire  frequency 
range  is  not  fully  utilized.  Note  also  that  the  sources  may  themselves  not  contain 
much  energy  at  low  frequencies.  Thus,  a  method  that  is  able  to  e.xploit  the  full  range 
of  available  frequencies  is  more  attractive,  and  the  next  subsection  presents  one  such 
method. 

Source  Location  Estimation  Method 

Instead  of  obtaining  source  location  estimates  from  the  relative  path  differences,  this 
method  obtains  the  source  location  estimates  from  }  directly,  and  hence  avoid 

the  problem  of  phase  wrapping. 

Consider  the  sensor  pair  (i,i)  located  at  the  same  radial  distance  and  separated 
by  an  angle  of  180  degrees,  see  for  example  Fig.  7.3.  This  type  of  sensor  pair  can  be 
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Figure  7.3:  Sensor  Pair 


found  in  common  array  configurations  such  as  ULAs  (uniform  linear  arrays),  circular 
arrays  and  cross  arrays.  Then  their  path  differences  (relative  to  the  reference  sensor 
at  the  origin)  are 


Of; 


dik  =  aisin(6k  +  <^,)  -  ■^cos'^(0k  +  0i). 

Ink 

=  -o^isiniOk  +  4>i)  -  r-^cos^iOk  +  0,). 

2/tfc 


(7.8) 

(7.9) 


where  Rk  denotes  the  range  of  the  fc-th  source  relative  to  the  reference  sensor,  Ok  is 
the  DOA  angle  of  the  fc-th  source  relative  to  the  reference  sensor  (measured  from  the 
vertical  axis),  a,-  is  the  radial  distance  of  the  i-th  sensor  from  the  reference  sensor  and 
0,  is  the  angle  position  of  the  i-th  sensor  relative  to  the  reference  sensor  (measured 
from  the  horizontal  axis). 


Both  the  path  differences  in  the  (i,i)  sensor  pair  have  the  same  curvature  term 
^cos^(0k  +  0t).  This  fact  is  used  in  decoupling  the  search  for  the  DO.A  angles  from 


the  search  for  the  source  range. 


Suppose  that  the  application  of  blind  identification  at  a  frequency  of  w  yields  the 
following  for  the  steering  elements  of  the  (/,/)  sensor  pair. 

=  exp{i^a,s:n(0,  +  0.)}  J 

We  ignore  the  gain  terms  since  they  appear  in  the  magnitudes  of  the  steering  elements 
<ind  can  be  normalized  to  unity  once  estimated.  Then  consider  the  problem  of  fitting 
the  steering  vector  for  the  {i,i)  sensor  pair,  say  c,(0,  R),  with  the  estimated  blind  id. 


estimate,  using 


where 


M9,R)=\\  c]i0.R)a,{0,.Rk) 


exp{-j^Qisin{e  +  0,)} 


c,{e,R)  =  -J  ^  ‘  exp{j-^cosHd + 0,)}  (7.12) 

[  cxp{;^a,szn(5  +  0.)}  ^^c2R  ^ 

Note  however  that  the  terms  involving  R  and  Rk  are  eliminated  in  the  above  fitting 

function  since  they  appear  as  scaling  phases  and  have  unit  magnitude.  The  function 

is  thus  dependent  only  on  the  angle  6  and  attains  a  maximum  value  at  the  angle 

9  =  9k.  For  L  sensor  pairs  we  define  the  fitting  function  as 


/(^)  =  II  il  cj(f?)a, 


where 


c.(^)  = 


^^p{-jjOiisin(0  +  <pi)} 


(7.13) 


(7.14) 


and  a,  denotes  the  blind  identification  estimate  of  a,{9k,Rk)-  The  DOA  estimate  for 
the  A:-th  source  is  then  found  by 


9k  =  arg  |max  f{0) 


(7.15) 


Remark;  The  Ar-th  source  range  can  be  estimated  once  the  A;-th  DOA  is  estimated 
by  defining  a  fitting  function  similar  to  the  one  above. 


104 


The  k-th  source  range  can  also  be  esiimated  independently  of  the  DOA  angle 
estimate  through  the  use  of  a  cross-array  structure,  see  Fig.  7.4  for  an  example  of 
a  cross-array.  In  this  case  the  search  for  the  range  can  be  conducted  in  parallel 


i^Sciuor 

Figure  7.4;  An  Example  of  a  Cross  Array 


to  the  search  for  the  DOA  angle.  Note  that  the  cross-array  can  be  viewed  as  con¬ 
sisting  of  two  sensor  pairs,  which  we  shall  denote  as  (?,  i)  and  with  locations 

({a,-,<^<},  {ai,0j-7r})  and  {{ai,<i>i  +  ir/2},  {a„<;i>,  -7r/2})  (in  polar  form),  respectively. 
Thus  the  path  differences  for  the  (i,i)  sensor  pair  is  as  found  in  Eqs.  (7.8),  (7.9)  while 
the  path  differences  for  the  (i*,i*)  sensor  pair  are  found  to  be 


d,.*  =  aiCos{6k  +  <i>i)-  TT^sin'^ldk  + 

Zlik 


(7.16) 


d,-,fc  =  -aiCOs(Qk  +  <i>i) 


(7.17) 
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(7.18) 


Writing  the  steering  eiennents  of  the  sensors  within  tlie  cross-array 

exp{-j^a{Sin{0k  +  Q,)}exp{j^j^cos^{0k  + 

a(5  R)=  +  +  (7  18) 

*'  exp{-j^QiCOs{9k  +  0i))eip{j'f-^sin^{9k  +  <P,)} 

exp{j^aiCOs{9k  +  0,)}exp{j^-^stn^i9k  +  0i)} 

Define  a  2-by-l  vector  hi{9k,Rk)  from  aj(0jt,/?fc)  by  the  following  procedure;  The 
first  element  in  is  set  equal  to  the  product  of  the  first  and  third  elements 

of  Sii{9k,Rk)  and  the  second  element  in  b,(flfc,  Rk)  is  set  equal  to  the  product  of  the 
conjugates  of  the  second  and  fourth  elements  of  a,(^jt, /?*)•  Writing  out  this  vector 
explicitly, 


exp{—j—ai{sin(0k  +  <Pi)  +  cos(9k  +  0,)))  (7.19) 

c 


Thus  the  range  has  been  decoupled  from  the  DOA  angle  and  using  a  similar  fitting 
function  as  previously,  we  estimate  the  range  for  the  A:-th  source  as 


where 


Rk  = 


g{R)  =11  d!(fl)b, 


and  b,  denotes  the  blind  identified  estimate  of  h,(0kk  Rk)  while 


MR)  =  f 

L  exp{-jj^ 


(7.20) 


(7.21) 


(7.22) 


Thus  the  estimation  for  the  DOA  angle  and  range  of  the  A:-th  source  can  be  performed 


in  parallel  since  they  are  decoupled  from  each  other. 


7.3.3  Source  Location  Estimation  using  a  Range  of  Fre¬ 
quencies 

At  this  stage  of  the  proposed  approach,  a  clustering  technique  is  applied  to  the  mul¬ 
tiple  estimates  of  source  locations  and  gains  found  for  a  given  range  of  frequencies. 
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The  final  source  locations  and  gains  estimates  are  computed  as  the  averages  of  the 
estimates  found  within  the  clusters. 

There  is  a  question  of  why  processing  over  a  range  of  frequencies  would  be  ad¬ 
vantageous  compared  to  that  of  just  using  a  single  bin.  The  obvious  answer  would 
be  the  fact  that  more  information  is  available  through  the  use  of  multiple  frequency 
bins  as  opposed  to  using  just  a  single  bin.  Furthermore,  clustering  of  the  multiple 
source  estimates  found  over  a  range  of  frequencies  is  a  natural  way  to  consolidate 
the  information  contained  in  the  given  range  of  frequencies.  The  reasons  for  this  are 
twofold. 

First,  at  any  particular  frequency,  the  sources  may  not  be  all  present.  Thus,  in 
order  to  estimate  all  the  source  locations  and  gains,  it  necessary  to  consider  more  than 
just  a  single  frequency  bin.  Suppose  that  at  a  particular  frequency  bin,  K  sources  are 
present.  Then  only  the  source  locations  of  the  K  sources  present  can  be  estimated. 
At  the  next  frequency  bin,  say  M  sources  are  present  and  their  source  locations 
estimated.  The  sources  at  this  frequency  may  or  may  not  be  found  in  the  previous 
frequency.  Note  however  that  the  sources  that  are  present  at  both  frequencies  would 
give  source  location  estimates  that  are  close  to  each  other  (ideally,  the  same  source 
locations),  while  the  location  estimates  of  the  sources  that  are  absent  at  the  previous 
frequency  would  be  further  away  (disjoint  from  the  other  sources).  By  sweeping 
through  a  range  of  frequencies  and  clustering  the  resulting  source  location  estimates, 
then  all  the  source  locations  may  be  found  from  the  clusters. 

The  second  reason  for  clustering  has  to  do  with  the  effect  of  observation  noise  on 
the  estimates  and  the  varying  relative  source  signal  powers  over  a  frequency  range. 
The  quality  of  data  at  different  frequencies  may  vary  significantly  over  the  range  of 
frequencies  on  which  the  blind  identification  algorithm  is  applied.  This  is  because  that 
the  source  signals  may  have  different  power  spectral  characteristics,  thus  varying  the 
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relative  energies  of  the  source  signals  at  clifTereni  frequencies.  Note  however  that  at 
the  frequencies  where  the  signal  powers  are  strong  (i.e..  high  SNRs),  the  corresponding 
source  location  estimates  found  through  blind  identification  should  be  near  to  the  true 
parameter  values.  The  frequencies  where  the  signal  energies  are  small  or  almost  non¬ 
existent  would  yield  parameter  estimates  that  are  erratically  distributed  (the  ‘‘outlier 
estimates”)  since  the  information  contained  within  those  frequencies  are  strongly 
dominated  by  the  spurious  observation  noise  power.  Note  that  the  noise  components 
at  different  frequencies  are  mutually  uncorrelated  and  therefore  it  is  not  likely  that  the 
“outlier  estimates”  would  map  to  the  same  locations  for  different  frequencies.  Thus 
the  application  of  blind  identification  over  a  range  of  frequencies  would  yield  clusters 
of  estimates  which  are  near  the  true  parameter  values  and  other  ‘outlier'  estimates 
which  are  erratically  distributed.  It  is  clear  then  that  any  post  processing  of  the 
multiple  source  location  estimates  would  neccessitate  the  omission  of  the  “outlier” 
estimates  which  fall  outside  these  clusters. 

We  present  one  clustering  technique  in  the  following  which  is  simple  and  have 
found  to  work  well  in  simulations.  Note  however,  that  there  are  various  other  clus¬ 
tering  techniques  found  in  pattern  recognition,  see  [37],  which  may  be  used  at  this 
point. 

One  clustering  method  is  by  the  use  of  a  2-D  histogram  of  the  multiple  estimates 
of  the  source  locations.  The  modes  (i.e.,  the  bins  within  the  histogram  where  the 
highest  count  of  source  location  estimates  fall  into)  of  the  histogram  is  selected  as  the 
center  of  the  clusters.  This  method  is  summarized  below. 

Clustering  via  2D  Histogram 

•  Collect  all  the  source  location  estimates  into  a  2-D  histogram.  The  number  of 
source  signals  can  be  found  by  determining  the  number  of  clusters  present  in 


108 


the  2-D  histogram  of  the  source  location  estimates.  Find  their  modes. 

•  For  a  specified  cluster  size  (A0  and  A/2)  of  the  modes,  pick  up  all  the  source 
location  estimates  falling  within  such  a  cluster  size  and  their  corresponding 
gains.  Compute  the  averages  of  these  estimates. 

This  averaging  of  the  estimates  falling  within  a  cluster  may  be  viewed  as  a  weighted 
averaging  of  all  the  parameter  estimates  found  in  the  range  of  frequencies  of  interest. 
That  is, 

for  the  A:-th  source  location 

=  (7.23) 

for  the  A;-th  source  gains 

where 

^(w;));  Source  location  estimate  found  at  frequency  wj. 
g(ufi):  Vector  of  gain  estimates  found  at  frequency  u./. 

.A/*:  fc-th  source  location  cluster. 

/<*:  Indicator  function  defined  as  follows 

1,  if  ((?(u;(),/2(u;j))  € 

0  otherwise. 

Here  Nk  =  53/  hk. 

Remark:  Other  weighting  strategies  may  be  employed  instead  of  the  indicator  func¬ 
tion  weights  given  above,  such  as  one  employing  variable  weights  which  give  greater 


(7.26) 
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weight  to  the  estimates  closer  to  the  cluster  centers  and  less  weight  to  estimates 
further  away. 

Note  that  the  post- processing  of  the  estimates  through  clustering  as  proposed 
above  is  fundamentally  different  from  the  coherent  approach  as  in  [75]  and  the  inco¬ 
herent  approach  as  in  [78].  Both  the  approaches  in  [75.  78]  applies  pre-processing  of 
the  observation  data  and  not  post-processing  of  the  parameter  estimates.  In  [75],  the 
observation  data  at  different  frequencies  is  focused  onto  a  single  frequency  while  in 
[78],  the  observation  data  is  averaged  across  different  frequencies. 

7.3.4  Source  Separation 

This  is  the  stage  where  the  final  source  location  and  channel  gains  estimates  is  used 
to  extract  the  source  signals  from  the  observation. 

The  source  signal  vector  s(t)  can  be  estimated  in  a  number  of  ways.  The  most 
straightforward  way  is  to  use  the  estimates  of  the  gains  and  source  locations,  {o,*,  0*,  /?*}, 
to  solve  for  S(a;),  i.e., 

S(u>)  =  A(w)*X(u;)  (7.27) 

An  alternative  approach  to  extract  the  different  source  signals  using  the  estimated 
source  locations  and  gains  would  be  to  use  spatial  filters  (beamforming  approach)  to 
extract  the  individual  signals,  see  [83]  for  the  various  beamforming  approaches  that 
may  be  applied. 

7.3.5  Summary  of  Proposed  Approach 

In  summary,  the  proposed  approach  proceeds  in  the  following  manner. 


Proposed  Approach 
1.  Apply  FFT  to  data  x(t). 
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2.  At  the  current  irequency,  estimate  the  number  of  sources  present  in  the  obser¬ 
vation.  This  can  be  accomplished  by  examining  tlie  singular  values  of  the  covarince 
matrix  or  through  such  criteria  as  found  in  [77,  26]. 

3.  Apply  Blind  Identification  to  X(u;).  The  locations  of  the  signal  sources  are 
estimated  from  the  estimated  matrix  A(u;)  and  the  channel  gains  are  estimated  from 
the  magnitudes  of  the  same  matrix.  Collect  these  estimates  of  gains  and  source 
locations. 

4.  Go  to  the  next  frequency  and  return  to  Step  2,  unless  the  maximum  desired 
frequency  has  been  exceeded.  If  so,  go  to  Step  5. 

5.  Apply  a  2-D  clustering  technique  onto  the  collection  of  source  locations.  Such 
clustering  methods  can  be  a  simple  histogram  or  the  methods  as  found  in  pattern 
recognition,  [37].  The  number  of  sources  is  estimated  from  the  number  of  clusters. 
Using  these  clusters,  estimate  the  averages  of  the  various  source  locations  and  asso¬ 
ciated  gains. 

6.  Using  estimates  of  {aik,6k,  Rk),  extract  {5i(u?)}  from  X{u;),  i.e,  S{t4;)  = 
A{u;)*X(u;). 

7.  Inverse  FFT  to  obtain  s(t). 

Discussions 

Unknown  Sensor  Gains:  As  is  clear  from  the  previous  section,  the  proposed  ap¬ 
proach  is  able  to  estimate  the  unknown  sensor  gains. 

Unknown  number  of  signals:  The  number  of  sources  is  estimated  after  clustering 
is  performed  on  the  source  location  estimates.  Thus  the  total  number  of  sources 
present  over  the  whole  range  of  frequency  bins  is  determined.  This  is  as  opposed  to 
the  estimation  of  number  of  sources  at  any  particular  frequency  bin  which  may  be 
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less  than  the  total  number  of  sources.  As  an  example,  consider  the  case  where  there 
are  four  sources,  two  of  which  have  overlapping  spectra  disjoint  from  the  other  two. 
Then  any  estimation  of  the  number  of  sources  at  any  particular  frequency  would  yield 
less  them  four  sources. 

Unknown  Source  Spectral  Characteristics:  The  proposed  approach  is  not  de¬ 
pendent  on  the  spectral  characteristics  of  the  sources.  Thus,  the  sources  can  be 
overlapping  or  non-overlapping  in  their  power  spectra. 

Unknown  Combinations  of  Narrowband  and  Wideband  Sources:  The  sources 
may  be  a  combination  of  narrowband  and  wideband  signals.  Our  approach  is  able 
to  still  cluster  the  narrowband  sources  even  though  the  number  of  points  within  the 
narrowband  source  cluster  is  relatively  less  than  the  number  of  points  within  the 
wideband  sources.  This  is  achieved  through  the  use  of  time  segmentation  of  the 
observation  data,  i.e.,  by  exploiting  the  fact  that  the  number  of  points  within  the 
narrowband  cluster  may  be  increased  by  using  multiple  observation  intervals.  Note 
that  the  ’outlier’  estimates  would  tend  to  distribute  themselves  erratically  with  dif¬ 
ferent  observation  intervals  while  the  ’good’  estimates  would  tend  to  stay  near  the 
true  estimates. 

Unknown  Combinations  of  Near-Field  and  Far-Field  Sources:  Source  loca¬ 
tion  estimation  is  performed  for  all  sources,  regardless  of  whether  they  are  far-field  or 
near-held.  Thus  the  fact  that  the  sources  may  consists  of  both  far-field  and  near-field 
sources  is  immaterial  to  the  proposed  approach. 
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7.4  Simulations 

7.4.1  Simulation  1:  Three  Near  Field  Sources 

The  array  is  a  cross  array  consisting  of  five  sensors  located  at  (0,0)  (reference  sensor), 
(.5,0),  (-.5,0),  (0,.5),  (0,-.5)  meters.  Three  source  signals,  one  of  interfering  white 
noise  with  power  equal  to  2.0,  while  the  others  are  a  male  speech  sequence  and  a 
female  speech  sequence,  both  of  unit  power.  The  interfering  white  noise  signal  is 
located  at  a  range  of  5  meters  and  angle  20  degrees,  the  male  speech  signal  is  located 
at  a  distance  of  1.5  meters  and  angle  -5  degrees  while  the  female  speech  sequence  is 
located  at  a  distance  of  1.5  meters  and  angle  20  degrees.  Additive  observation  noise 
is  taken  to  be  AWGN  (additive  white  gaussian  noise)  with  variance  equal  to  (.1)^. 
Speed  of  propagation  is  taken  to  be  340  m/s. 

The  channel  attenuations  in  the  model  is  assumed  to  be  due  only  to  geometrical 
attenuation  of  the  signal  sources  (from  the  spherical  spreading  of  the  wavefront)  and 
thus  the  gains  a,fc  are  inversely  proportional  to  the  distances  i?,*,  i.e.,  the  distance 
between  the  i-th  sensor  and  the  k-ih  source.  The  gains  are  computed  to  be  = 
R\klRik<t  where  the  gains  for  the  first  sensor  are  normalized  to  unity.  Then  the  gains 
are 
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16000  samples  of  observations  from  the  array  is  taKen  at  a  sampling  rate  of  8000 
Hz. 

These  samples  are  then  chopped  into  30  blocks  of  non-overlapping  intervals  of  512 
samples  each  and  512  point  FFT  is  applied  to  these  intervals.  Wideband  blind  identi¬ 
fication  is  then  applied  on  the  data  for  frequency  bins  ranging  from  50  Hz  to  1 .6  kHz. 
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Figure  7.5:  Simulation  1:  Typical  observations  from  sensors  no.  1  and  2. 

Note  that  this  selection  of  frequency  bins  does  not  satify  the  constraint  for  avoiding 
phase  wrapping  and  therefore  the  closed  form  least  squares  based  methods  cannot  be 
applied.  The  source  locations  are  estimated  using  the  Approximate  Method.  Figs. 
7.6  show  the  2-D  histogram  plots  of  the  source  location  estimates. 

It  is  clear  from  the  2-D  histogram  that  there  are  three  definite  clusters  and  there¬ 
fore  there  are  three  souce  signals  present. 

Picking  the  size  of  the  clusters  to  be  A0  =  5  deg  and  Ai?  =2  meters,  i.e.,  all 
source  location  estimates  falling  within  this  distance  from  the  modes  of  the  clusters 
are  collected  and  averaged,  we  find  that  the  final  estimated  DOA  angles,  ranges  and 
gains  are:  (19.79  deg,  4.21  meters),  (-4.48  deg,  1.65  meters)  and  (19.15  deg.,  1.83 
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figure  7,6:  Simulation  1:  2-D  Histogram  of  the  Source  Location  Estimates, 


meters),  respectively,  and  their  corresponding  gains  are 
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(7,29) 


Using  these  averaged  gains  and  source  location  estimates,  the  source  signals  are  ex¬ 
tracted  from  the  observation  data.  The  true  and  separated  speech  sequences  is  shown 
in  Figs,  7.7-7. 8.  Computation  of  the  correlation  between  the  estimated  and  the  true 
speech  sequences  shows  the  correlation  to  be  99  percent  and  95  percent  respectively. 
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Figure  7.7:  Simulation  1:  Plots  of  Estimate  of  Speech  no.  1. 

7.4.2  Simulation  2:  Near-Field  (Wideband)  Source  and  Far- 
Field  (Narrowband)  Source 

The  shape  of  the  array  is  the  same  as  in  Simulation  1,  except  that  the  radial  distances 
of  the  sensors  from  the  reference  is  set  to  0.25  meters.  The  wideband  source  is  a  speech 
source  located  at  range  of  1.5  meters  and  angle  -45  degrees.  The  narrowband  source 
is  synthetically  generated  to  contain  sinusoidal  lines  at  the  frequencies  of  55.  60  and 
65  Hz.  It  is  situated  at  a  range  of  20  meters  and  angle  of  5  degrees.  Thus  we  have 
a  combination  of  wideband  near-field  source  and  narrowband  far-field  source  in  the 
observations.  Both  sources  are  at  unity  power.  Additive  observation  noise  is  taken 
to  be  AWGN  with  variance  equal  to  (.1)^.  Speed  of  propagation  is  taken  to  be  340 
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Figure  7.8:  Simulation  1:  Plots  of  Estimate  of  Speech  no.  2. 

m/s. 

The  channel  attenuations  in  the  model  is  assumed  to  be  due  only  to  geometrical 
attenuation  of  the  signal  sources  (from  the  spherical  spreading  of  the  wavefront)  and 
thus 
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As  done  previously,  16000  samples  are  collected  and  30  sub-intervals  of  512  points 
each  are  obtained.  The  range  of  frequencies  over  which  the  proposed  approach  is 
applied  is  from  15  Hz  to  234  Hz,  giving  a  total  of  15  frequency  bins.  The  proposed 
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Figure  7.9:  Simulation  2:  Typical  observations  from  sensors  no.  1  and  2. 

algorithm  is  applied  to  the  FFT  of  the  data.  Then,  another  30  sub-intervals  of  512 
points  each  is  obtained  from  the  same  observation  set  by  shifting  the  sub-intervals 
window  by  200  points  and  the  proposed  approach  is  again  applied.  Next,  the  same 
procedure  is  done  for  a  sub-interval  window  shift  of  400  points,  and  also  for  a  shift 
of  600  points.  Finally,  the  collection  of  source  location  estimates  is  clustered  using  a 
2D  histogram,  shown  in  Fig.  7.10.  Two  clusters  are  evident  from  the  histogram. 

Picking  the  size  of  the  clusters  to  be  A©  =  5  deg  and  Ai?  =2  meters,  i.e.,  all 
source  location  estimates  falling  within  this  distance  from  the  modes  of  the  clusters 
are  collected  and  averaged,  we  find  that  the  final  estimated  DOA  angles,  ranges  and 
gains  are:  (-44.6  deg,  1.43  meters)  and  (5.23  deg.,  25  mete.  ),  respectively,  and  their 
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Figure  7.10:  Simulation  2:  2-D  Histogreim  of  the  Source  Location  Estimates. 


corresponding  gains  are 
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Using  these  averaged  gains  and  source  location  estimates,  the  source  signals  are  ex¬ 


tracted  from  the  observation  data.  The  true  and  separated  speech  and  narrowband 


sources  are  shown  in  Figs.  7.11-7.12. 


Figure  7.11:  Simulation  2:  Plots  of  Estimate  of  Speech  Source, 

7.4.3  Simulation  3:  Far-Field  Sources 

It  is  assumed  that  the  number  of  sources  is  three  and  the  number  of  sensors  (i.e., 
microphones)  is  four.  The  sources  are  two  male  speech  signals  and  a  female  speech 
signal,  respectively  with  no  multipaths.  The  sources  are  located  in  the  far  field  with 
DO  A  ajigles  of  -40,  +10  and  +50  degrees  respectively.  The  four  sensors  are  located 
at  (0,0),  (1,0),  (1, 1),  (0, 1)  meters  respectively.  The  speed  of  propagation  is  taken  as 
340  meters- per-second. 
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Figure  7.12:  Simulation  2:  Plots  of  Estimate  of  Narrowband  Source. 


The  gains  are  arbitrarily  selected  as  follows: 
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(7.32) 


16,000  time  samples  of  the  microphone  array  are  taken  at  a  sampling  rate  of  8000 
samples/s.  The  signal  sources  are  all  of  equal  power,  with  the  total  energy  of  a  signal 
source  sununed  over  the  16,000  samples  set  to  unity.  The  separation  of  the  speech 
signals  from  the  array  observations  is  desired.  The  time  samples  are  divided  into  30 
sub-intervals  of  512  samples  each  and  512-pt  FFT  taken  over  the  3U  sub-intervals. 
Wideband  blind  identification  is  then  applied  across  the  frequency  spectrum  ranging 
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from  93.75  Hz  to  703.125  Hz,  giving  a  total  of  40  frequency  bins.  DO  A  estimates  are 
found.  The  incoherent  clustering  algorithm  is  applied  to  obtain  the  final  DOA  and 
gain  estimates.  Using  these  estimates,  individual  estimates  of  the  speech  sources  are 
made,  using  (7.27). 

Modes  of  DOAs  are  found  to  be  at;  -41  deg.,  +10  deg.,  and  +50  deg.  A  plot  of 
the  histogram  of  the  DOA  estimates  is  shown  in  Fig.  7.13. 


Figure  7.13:  Simulation  3:  Histogram  of  DOA  estimates 

Picking  the  size  of  the  clusters  to  be  A0  =  5  deg,  i.e.,  all  ,  we  find  that  the  final 
(incoherent)  estimated  DOA  angles  are: 

Est.  DOAs:-40.667  deg,,  +10.417  deg.,  +50.0  deg. 

Also,  using  the  sensor  gains  correponding  to  the  collected  DOA  estimates,  we  find 


that  the  hnal  gain  estimates  are: 
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1.06 

1.14 

1.794 

1.23 

1.572 

(7.33) 


Figs.  1  (a,b)  show  two  typical  observations  from  the  array.  It  can  be  seen  from  Figs. 
2(a,b),  Fig.  3(a,b),  Fig.  4(a,b)  that  the  speech  signals  are  succesfully  separated  from 
the  array  observations. 


True  speech  simal  no.  1 
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Figure  7.15:  Simulation  3:  Estimate  of  Speech  no.  2 

7.4.4  Simulation  4:  Non-stationary  Source  and  Speech  Source 

The  scenario  is  the  same  as  is  in  Simulation  2,  except  that  the  narrowband  source 
is  now  replaced  by  a  non-stationary  source.  The  non-stationary  source  is  generated 
by  scaling  a  gaussian  distributed  source  differently  at  each  time  instant,  thus  making 
the  source  have  different  ensemble  variances  at  different  times. 


The  proposed  approach  is  applied  over  a  frequency  range  of  50  Hz.  to  1.6  kHz. 
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Figure  7.16:  Simulation  3:  Estimate  opeech  no.  3 


True  Gains 


■  1.0 

1.0 

’  1.0 
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1.12 

1.01 

1.07 

.927 

.89 

1.0 

Estimated  Gains 

.89 

.984 

.89 

.99 

.943 

1.03 

1.12 

1.0 

1.11 

.992 

(7.34) 


The  true  source  locations  are:  (1.5  m,  -45  deg.)  and  (20  m,  5  deg.) 

The  estimated  source  locations  are:  (1.31  m,  -45.08  deg.)  and  (25  m,  4.393  deg.). 

From  fig.  7.17,  it  is  clear  that  the  speech  signal  is  succesfuly  separated  from  the 
non-stationary  signal. 
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Figure  7.17;  Simulation  4:  Typical  observation  and  estimated  speech  source. 

7.4.5  Simulation  5:  Non-stationary  Source  and  Two  Speech 
Sources 

The  scenario  is  the  same  as  is  in  Simulation  1,  except  that  the  white  noise  interference 
source  is  now  replaced  by  a  non-stationary  source.  The  non-stationary  source  is 
generated  by  switching  between  uniform  and  gaussian  distributions  every  250  time 
samples  where  their  variances  are  different. 

The  proposed  approach  is  applied  over  a  frequency  range  of  50  Hz.  to  1.6  kHz. 
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Figure  7.18:  Simulation  5:  Typical  observations  from  sensors  no.  1  and  2. 


True  Gains: 
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(7.35) 


Using  these  averaged  gains  and  source  location  estimates,  the  source  signals  are  ex¬ 


tracted  from  the  observation  data. 


Computation  of  the  correlation  between  the  estimated  and  the  true  speech  se¬ 
quences  shows  the  correlation  to  be  99  percent  and  94  percent  respectively. 
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Figure  7.19:  Simulation  5:  2-D  Histogram  of  source  location  estimates. 

7.5  Extension  of  Approach  to  Multipath  Cases 

The  proposed  approach  is  also  applicable  to  the  multipath  case,  specifically  the  spec¬ 
ular  multipath  scenario,  i.e.,  perfectly  coherent  multipaths  are  present.  In  this 
subsection,  we  will  discuss  the  application  of  wideband  blind  identification  to  this 
problem. 

For  ease  of  discussion,  assume  1  multipath  signal  and  n  other  signads. 

The  observation  at  the  i-th  sensor  is  written  as 


^.(0  =  ULi  <^ikSk{t  -  Tik)  +  -  6)  +  m- 


$ 
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Figure  7.20:  Simulation  5:  Plots  of  Estimate  of  Speech  no.  1. 

denotes  the  multipath  ga.ns  and  delays.  Thus  the  signal  si(0  has  a  direct 
and  a  multipath  component. 

In  the  frequency  domain, 

X(u;)  =  A(u;)S(u,)  +  vi^) 


where 


A(a;)  = 


1.0 


(7.36) 
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Figure  7.21:  Simulation  5:  Plots  of  Estimate  of  Speech  no.  2. 

Blind  identification  is  then  performed.  The  signal  with  multipath  can  be  distinguished 
from  the  other  signals  (with  no  multipaths)  by  examining  the  column  vectors  of  the 
matrix  A(u;).  The  magnitudes  of  elements  within  the  column  vector  associated  with 
the  signal  si{t)  would  vary  with  different  u>'s  while  those  with  no  multipaths  would 
have  constant  magnitudes.  Therefore,  the  estimates  of  for  t  =  2,...,m, 

k  =  2,...,n  can  be  found  directly  from  A.  Estimates  of  {aix,r,i}  and  may  be 

found  from  the  estimates  of  A(w)  at  several  different  frequencies. 

In  general  however,  the  use  of  blind  identification  raises  the  problem  of  scaling 
2md  ordering  indeterminacy  of  the  signals,  see  (70).  When  there  cure  several  signals 
with  multipaths  this  indeterminacy  is  undesirable  because  in  order  to  solve  for  the 
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direct  and  multipath  parameters,  estimates  of  A  must  be  made  at  several  frequencies. 
However  with  several  signals  with  multipaths,  their  corresponding  column  vectors 
in  A  may  be  scaled  and  ordered  differently  at  different  frequencies.  Since  signals 
with  no  multipath  can  be  easily  distinguished  from  signals  with  multipath  (as  noted 
previously)  it  is  assumed  that  A  contains  only  those  column  vectors  corresponding  to 
signals  with  multipaths.  Without  loss  of  generality,  we  consider  the  case  where  there 
are  two  signals  each  with  a  multipath.  Thus, 


(L0  +  ;9iie--'“'<«‘)  (1.0  + 

)  (a22e"-'‘*^*  +  /322e"-’“^” ) 


(7.37) 


To  eliminate  the  unknown  scaling  of  the  column  vectors  in  the  blind  identified  value 
of  (7.37),  a  ’normalized’  matrix  A(u;)  is  defined  by  dividing  the  column  vectors  of  the 
estimate  of  the  matrix  with  their  first  elements,  i.e.. 


1.0  1.0 

l+/3ue-/“hi  l+ftje  •'■'‘I* 


(7.38) 


Thus  the  unknown  scaling  of  the  column  vectors  in  A  at  different  frequencies  have 
been  eliminated.  .Next  the  problem  of  ordering  of  the  column  vectors  is  examined. 


Define 


/.(^)  - 


+  /?2.e--''-<»*) 


(7.39) 


When  the  separation  between  adjacent  bin  frequencies  (Aw)  are  small  enough  such 
that  I  Awr^fe  |<<  1.0,  |  Aw^ifc  1<<  1.0,  then 


/.(wo  +  Aw)  =  /i(wo)  -  jAwgi(wo) 


(7.40) 


where, 


(fl2<r2.e-J‘^> 

1  + 


(7.41) 


Similarly  at  the  bin  frequency  u;o  +  Aw,  we  have 


/,(wo  -  Aw)  =  /,{wo)  +  j  Awfif,(wo )  ( ”-42) 

Now  consider  the  ratio  of  differences: 

^  /i(u,  +  Au.)  -  hM  J, 

//(wo  +  Aw)  —  /t(wo) 

it  can  be  seen  from  (7.40)  and  (7.42)  that  generally,  piki  is  equal  to  -1.0  for  z  =  fc  =  I 
and  not  equal  to  -1.0,  otherwise.  Thus,  given  the  matrices  A{wo),  A(wo  +  Aw), 
A(wo  —  Aw)  (with  possibly  different  and  unknown  permutations  of  their  column 
vectors),  the  ratio  in  (7.43)  can  be  used  to  determine  their  correct  ordering. 

Thus  the  re-scaling  and  re-ordering  of  the  column  vectors  may  be  accomplished 
as  follows: 

Scaling:  Define  the  ’normalized’  matrix  A(w)  as  in  (7.38). 

Ordering:  Given  the  matrices  A(wo),  A(wo  +  Aw),  A(wo  -  Aw). 

1.  Suppose  the  correct  column  vectors  in  A(wo  +  Aw),  A(wo  —  Aw)  are  to  be 
associated  with  the  z-th  column  vector  in  A(wo).  Denote  the  second  element  in  the 
z-th  column  vector  of  A(wo)  as  o. 

2.  Denote  the  elements  in  the  second  row  of  A(wo+ Aw),  A(wo-'Aw)  as 

and  •••}  respectively.  Here  the  term  pk  denotes  the  second  element  of  the  k-th. 

column  vector  in  A(wo+ Aw)  and  similarly  for  Uk.  Define  the  differences  Spk  =  Pk—oi^ 
and  Suk  =  Vk  ~  ct. 

3.  Compute  the  ratios:  |^.  The  ratio  closest  to  the  value  of  -1.0  would  signify 
the  correct  ordering  of  the  column  vectors,  i.e.,  the  {j,  k)  index  of  the  ratio  denotes 
that  the  7-th  column  vector  of  A(wo-j- Aw)  and  the  k~th  column  vector  of  A{wo  — Aw) 
are  associated  with  the  z-th  column  vector  of  A(wo). 
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7.6  Extension  of  Approach  to  3-D  Scenario 

The  proposed  approach  may  also  be  extended  to  the  case  where  3  parameters  deter¬ 
mine  the  source  location  (eg.,  range,  azimuth  angle  and  elevation  angle).  Conside  fig. 
7.22.  The  path  differences  (relative  to  the  reference  sensor  which  is  assumed  to  be  at 

k'th  lource 


Figure  7.22:  3-Dimensional  Scenario  for  Sources  and  Sensors 
the  origin)  can  then  be  shown  to  be 

A*  =  4  -  -  (dL?)  (T.44) 

where  =  a,{sin7fc sin^< sin (0fc  +  ^,)},  is  the  far-field  component  of  the  time 
delay. 

Similarly,  to  the  previous  discussion  on  the  sensor  pairs  in  the  2-D  case,  the  estima¬ 
tion  of  the  azimuth  and  elevation  angles  may  be  decoupled  from  the  estimation  of  the 
range  through  the  use  of  a  specific  sensor  pair  configuration.  Specifically,  suppose  that 
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the  z’-th  sensor  aissociated  with  the  i-th  is  located  at  the  location  (a,,  — tt  -f  <g,  tt  — 

c 

Then  its  path  difference  for  the  t-th  source  is  written  as 

=  (T-IS) 

where  <i,4  is  as  defined  previously  for  the  i-th  sensor.  Thus,  by  defining  a  fitting 
function  similar  to  the  one  in  the  2-D  case,  the  search  for  the  azimuth  and  elevation 
angles  {6k',lk)  is  decoupled  from  the  range  search. 

Remarks: 

(i)  The  search  for  the  source  range  may  be  performed  once  the  azimuth  and 
elevation  angles  are  estimated. 

(ii)  On  the  other  hand,  it  may  be  possible  to  configure  the  sensor  array  structure 
analogously  to  the  cross-array  strucutre  of  the  2-D  case,  such  that  the  source  range 
can  be  estimated  independently  of  the  azimuth  and  elevation  angle  estimation. 


8 


CONCLUSIONS 


We  conclude  the  with  the  following  summary  of  the  results  of  this  investigation. 

•  An  analysis  of  ESPRIT  under  random  model  errors  is  conducted  which  provides 
insight  into  its  performance  and  robustness  against  such  random  model  errors 
as  unknown  sensor  gains  and  phases,  unknown  sensor  locations  and  unknown 
rotations  of  sub-arrays  within  the  ESPRIT  array.  In  particular,  it  is  shown  that 
in  the  case  of  Uniform  Linear  Arrays  (ULA),  the  bulk  of  the  Mean- Square- Error 
(MSE)  of  ESPRIT  is  due  to  the  sensor  phase  errors  as  opposed  to  the  sensor 
gain  errors. 

•  An  approach  based  on  a  signal  subspace  constraint  is  proposed  for  array  sig¬ 
nal  processing  under  model  errors.  This  approach  yields  asymptotically  unbi¬ 
ased  estimates  of  the  sensor  gains  and  phases  when  the  DOA  angles  are  known 
a-priori.  When  both  sensor  gains/phases  and  DOA  angles  are  unknown,  an 
iterative  version  of  the  approach  is  given. 

•  A  novel  approach  using  blind  identification  is  developed  for  the  source  estima¬ 
tion  problem  under  model  errors  and  uncertainties.  By  using  clustering  tech¬ 
niques  in  conjunction  with  blind  identification,  it  is  shown  that  the  proposed 
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approach  is  robust  against  the  silmutaneous  presence  of  such  uncertainties  such 
as  unknown  sensor  or  channel  gains,  unknown  combinations  of  near-field  and 
far-field  sources,  unknown  combinations  of  wideband  and  narrowband  sources, 
unknown  source  spectral  charateristics  and  unknown  number  of  sources.  As 
such,  it  shows  promising  potential  in  such  application  areas  as  speech  separa¬ 
tion/processing  and  mobile  communications. 
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